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ABSTRACT 


The Harris expansion method is applied to the elastic s-wave 
scattering of low energy electrons from hydrogen atoms and singly 
ionized helium atoms. The trial wave functions are Hylleraas 
functions of 22, 34, 50 and 70 parameters. It is found that rea- 
sonably accurate values of e-H phase shifts can be calculated but 
that e-He™ phase shifts are substantially less reliable. It is 
Shown that the Harris method gives an accurate depiction of the 
location, but not the width, of the scattering resonances. Singlet 
and triplet s-wave phase shifts for e-H and e-He” scattering are 
compared with the results of other calculations and H and He § 
state energy levels including the auto-ionizing levels are pre- 
sented and compared with other calculations and with experiment. 
It 1S tentatively concluded that the Harris method does not work 


on systems whose long range potential is of the Coulomb type. 





TABLE OF CONTENTS 


I. ie enn le ee ee ee ee ee ee ee 7 
tis EXPANSION TECHNIQUES IN SCATTERING THEORY----------- 11 
eee AL WAVES “ANePHASE SHIFTS<---------------=-- 11 

B. THE KOHN VARIATIONAL METHOD--------------------- 7 

C. THE HARRIS EXPANSION METHOD--------------------- 18 

1. Formal Solution----------------------------- 18 

2. Application to Atomic Scattering------------ p2 

3. Summary of the Harris Method--~-------------- 33 

III. RESULTS AND DISCUSSION ------9-%- once r rr rrr rrr eter ec rn- 35 
A. ELECTRON-HYDROGEN SCATTERING ---c-cc coco coerce 36 

l. Phase Shifts -------------------------------- 36 

2. Scattering Resonances ----------------------- 47 

2 Ce Sic See |S a ee ee 55 

B. ELECTRON-HELIUM ION SCATTERING ------------------ Dy) 
re 0 a 7 

2. Scattering Resonances ----------------------- 67 

Ce ey Ie veyo 2 ----=-- = — ay 

LEG CO GUSTO Se: = = Snes ae eae eo 7h 
APPENDIX A Formulae for the Various Matrix Elements ------ 76 
ar eretx B Imibearal Recursion Relations--------------==-<-- 82 
ewe € Kayliecigh-Ritz Variational Method -------------- 85 
APPENDIX D Discussion of Numerical Methods---------------- By 
7 ee eee i ee so So ee So ee Set 105 
REFERENCES ---~-------------------------------------------- 168 
INITIAL DISTRIBUTION LIST-----------~----------------------- 171 
FORM DD 1473 ------ nnn nn enn nr rer rn ener rrr rrr ncn 172 





ee. 


Lie 


IV. 


PrSt-On FABEES 


Beets Set Size for Various M=-==<..+~---<-.-. 2265-65 -s 
Singlet and Triplet Phase Shifts for e-H Elastic 
Scatter ing---------------------- -------- +--+ -- +--+ 
H Auto-ionizing Levels Corresponding to e-H 
Scattering Resonances -------------------+------------- 
Singlet and Triplet Phase Shifts for e-He” 

Elastic Scatter ing---------------------------------- 
He Auto-ionizing Levels Corresponding to Senet 


Scattering ResonancCeS - ---<-- <4 -- 2 ee en een eee eee 





aye 


aS; 


Loli Ons: TGURES 


Integration Region for A(\V,A ,2)--------------------------- 26 
Interchanging the Ty and I 5 eae erat 1G ia le 31 
Eigenvalue Trajectories for Singlet e-H System. N = 70----38 


Typical Evolution of One Eigenvalue Trajectory as N is 


Incr CAS C= wm er er en rr rn nn ee enn rn ee re een eee cne- 40 
Singlet Phase Shift for e-H as a Function of @. k = 0.4---42 
Singlet Phase Shift for e-H as a Function of @. k = 0.8---44 
Evaluation of e-H Singlet Scattering Length. N = 70------- L6 
Eigenvalue Trajectories for Triplet e-H System. N = 70----48 


Triplet Phase Shift for e-H as a Function of @. k = 0.4---49 


Triplet Phase Shift for e-H as a Function of @. k = 0.8---50 
Calculated Phase Shifts for e-I] Scattering---------------- oo 
Singlet Energy Levels for pe = eee ee ee 58 
Triplet Energy Levels for H . N = 50---------------------- 59 


E1lgenvalue Trajectories for Singlet e-He System. N = 70--61 
Singlet Phase Shift for e-He’ as a Function of @. k = 0.8-62 
Singlet Phase Shift for e-He as a Function of @. k = 1.4-63 
Eigenvalue Trajectories for Triplet e-He System. N = 50--64 
Triplet Phase Shift for e-He as a Function of @. k = 0.8-65 


: : ~ a : 
Triplet Phase Shift for e-He as a Function of @. k = 1.4-66 


Calculated Phase Shifts for ecltes Scattering-c--e--ee-- eH 69 
Singlet Energy Levels for He. N = 70---------------------- Ge 
Triplet Energy Levels for He. N = 50---------------------- 73 
Example of Method of Choosing h; (x) Ce a ere 98 





ACKNOWLEDGEMENT 


It 1S a pleasure to acknowledge the continued support and 
encouragment of Professor Otto Heinz and also that of Professor 
Robert Armstead der “hose direction this work has proceeded. 

Dr. D.E. Harrison, Jr., provided several helpful suggestions re- 
garding techniques of numerical integration. Professor A.H. Stroud 
of Texas A & M University Kindly provided the computer program by 
which the sample points and weight factors for the Gauss-Legendre 
numerical integrations were generated. The numerical calculations 
were performed on the IBM 360/67 computer of the Naval Postgraduate 
School Computer Facility, whose personnel were most helpful on 
several occasions. Last but not least, the author wishes to grate- 
fully acknowledge the continued moral support of his wife, who 


perservered through some very trying times. 





I. INTRODUCTION 


The quantum theory of scattering has had a long and fruitful 
history, beginning with the very origins of quantum theory. Its 
extenSive application to atomic collisions is indicated by the 
large number of texts, monographs and reviews of that subject pre- 
sently in print LiJ. More recently the wide availability of high 
speed electronic computers has given additional impetus to the 
Subject. New methods have been devised and older ones refurbished 
to take advantage of the computational power now at the disposal 
of virtually every investigator in the field. Parallel to the 
development of computer technology has been the application of 
powerful variational techniques to the scattering problem. 

Variational methods have held a distinguished place in mathe- 
matical physics since the days of the Bernoullis. Perhaps no 
technique has had wider application and thus it 1S no surprise that 
it has proven a powerful tool in scattering theory as well. Its 
first application to atomic scattering by Hulthen in 1944 (245 
followed shortly by Kohn (3] and Lippman and Schwinger [4] has led 
to extenSive applications to all areas of scattering theory (5). 
One difficulty that arose with the Kohn and related methods how- 
ever, was the appearance of spurious singularities in the cal- 
culated phase shifts. While not fatal to the method they were 
annoying and some considerable effort has been expended in attempts 
to explain these singularities by Schwartz (6] and others Ceol 


By working to exploit rather than eliminate the singularities, 








Harris [10] developed an expansion technique that has yielded very 
reliable results to date. While this scheme involves some compu- 
tational simplifications over the Kohn method it too is not without 
disadvantages, which are dealt with in some detail below. 

Formal solutions abound to the scattering problem, ranging 
from the Born approximation and partial wave analysis to the 
Green's function methods developed by Schwinger C11], but accurate 
solutions to "real'' problems beyond the realm of potential scat- 
tering are relatively few, owing principally to the complexity of 
the problem when one permits the scattering center to have 
structure. Even the simplest of "non-trivial" scattering problems, 
that of electrons or positrons incident elastically on ground 
state hydrogen atoms 1S already a "many body"’ problem and thus nee 
Subject to solution in closed form. 

It is the purpose of this thesis to investigate the application 
of the method developed by Harris to the problem of the elastic 
scattering of low energy electrons from hydrogen atoms and from 
Singly ionized helium. Both have been attempted several times 
previously by various other methods, most notably in the case of 
hydrogen by Schwartz, whose definitive calculation was published 
in 1961 [12] and more recently by Callaway and his collaborators 
in 1969 and 1970 asec. and by Adelman and Reinhardt [15] and 
Truhlar and Smith [16] in 1972. Because of the more comp lex nature 
of the interaction in the case of singly ionized helium and a 
paucity of experimental data for comparison these calculations 
have not been as extensively pursued. The earliest attempt was 
by Bransden and Dalgarno in 1953 [17] and the most recent by Burke 


and Taylor in 1966 [18]. 





With the exception of Schwartz' work most of the important 
calculations of these problems have used the algebraic close cou- 
pling approximation [19] wherein the trial wave function is made 
Mpeeatemeast oan part of a combination of the first few botind state 
orbitals of the target atom. Such a scheme has earned a deserved 
reputation for accuracy and rapid convergence, however from a 
philosophical point of view it suffers because its application 
requires some a priori knowledge of the detailed structure of the 
fam@et particle. On the other hand Schwarta, in his calculation, 
made use of an ever increasing set of those functions first used 
by Hylleraas in his famous calculation of the ground state energy 
of helium [20]. This method makes no approximations other than 
those required by the use of a finite rather than an infinite basis 
Space in which to work? 

It was in the spirit of Schwartz' work that the present project 
was undertaken. It was hoped that by utilizing the Hylleraas type 
functions in a calculation of electron-helium i10n scattering that 
a more accurate and complete compilation of s-wave phase shifts 
would be obtained. As will be seen below such hopes were not 


realized, however it is hoped that the reporting of such results 


cenniiitiemedtiiite ot 





I ohis Heenot Strictly true. Electron capture by the seattering 


center is usually neglected in calculations of this type, and is 
neglected also in the calculation reported in this thesis. To 
incorporate capture requires including terms in the potential in- 
volving the interaction of the electrons with the electromagnetic 
field in order to account for the creation of the necessary photons 
and iS properly a problem in quantum electrodynamics. The reported 
cross sections for capture in the system under consideration here 
are small compared to the elastic scattering cross sections and 
usually may be safely neglected. 








as were obtained will provide incentive for others to seek an 
explanation of the apparently anomalous results for helium ion 
scattering neces reported. 

In Section II a brief development of the partial wave method 
of scattering theory will be presented including the modifications 
required in the presence of the long range Coulomb force followed 
by a description of the Kohn method. Last will be a more detailed 
description of the Harris method and its application to elastic 
electron scattering by single-electron atoms and ions. In Section 
III the results of the numerical calculations are presented for 
both atomic hydrogen and singly 10nized helium along with a dis- 
cussion of these results and comparison with those of other workers. 
An important by-product of the Harris method is the calculation of 
the bound state energy levels of the corresponding two-electron 
System. Results of these calculations for H and He are also 
presented in Section III. Most of the detailed calculations as 
well as discussion of the numerical methods employed are deferred 


to the appendices. 
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II. EXPANSTON TECHNIQUES IN SCATTERING THEORY 


In this section will be presented a brief outline of the method 
of partial wave analysis [21] for the case of elastic atomic scat- 
tering. Part B will briefly present the Kohn method and in part C 
the Harris method will be developed in some detail along with its 
application to elastic He" scattering of electrons. The system of 
units used throughout will be those in which m (electron mass), 
kR (Planck's constant) and e (electronic charge) all equal unity. 
In such a system the unit of energy is equal to twice the ground 
state energy of hydrogen (1i.e., 27.2 eV), and momenta are written 
in terms of the wave number, k, such that E = aioe It will be 
assumed that the nuclear mass, be it hydrogen or helium, is infi- 


nitely large compared to that of the electrons. 


A. PARTIAL WAVES AND PHASE SHIFTS 

Let a beam of particles of momentum k fall upon a spherically 
symmetric potential V(r). Taking the scattering center as the 
coordinate origin and the direction of k as the angular axis the 
wave function of the outgoing beam far from the coordinate center 
will be proportional to the sum of the wave function of the inci- 
dent beam, represented by the plane wave exp(ik-r) and a modified 
spherical wave, f(@)exp(ikr)/r, representing the scattered portion 
of the incident beam. For proof that solutions to the Schroedinger 


equation of the form 


+ (onenes (1) 
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@east, the reader is directed to reference 21 (p. 220 ff). The 


differential scattering cross-section is defined to be 


Number of particles per unit time scattered into 
the solid anole element dn atin 


Oo = 
wid Total lnerdentr lux (2) 








where the particle flux is defined to be 


ye 4 * 
ase Y Vaca), (3) 





= — = A D, 
sae O_ % 16 
: 0 = = eee 
(9) dn 3 = da (4) 
O "oO 


—_—? , A - 
where ds 1S the element of surface area at the detector and sris a 


Manet Vector in the direction of the scattered flux. fo and ve 
represent respectively the flux arising from the first and second 
terms of (1), respectively. Because of the scalar product in (4) 
only the radial part of the scattered wave function contributes 

to the scattered flux. Evaluating es and wae from (3) and applying 


them to (4) leads immediately to 
= | |< 
O(O)dn = |£(O)} da . (5) 


It remains then to exhibit a means of calculating £(9) to 
complete the connection between theory and experiment. Such a 
means is the method of partial wave analysis. The most general 
solution to the Schroedinger equation in the case of axial sym- 
metry 1S 


wy (x) 
~~ y°(0) (6) 





er) =e 
£ 


ZZ 





| O , ’ 
where Y¥, (8) 1s the spherical harmonic of order zero. The index &£ 
represents the angular momentum quantum number. Each radial function 


satisfies an equation of the form 





2 Use) 
ia) 2450 £(L+1) k £ a 
:- Be oF (x? Se )+ ves ee | eames oe 


If V(r) goes to zero fast enough at large r there will be a region 
outside of which its effects can be neglected. Since uy (x) isa 


function of r only, equation (7) can be rewritten 


Hoa eee Jaye i ao 


and since the last two terms tend to zero for large r it is rea- 

sonable to assume that any solution u tends to the form A*SsSin(kr+t€) 
: aE) sie 

for large r. If this be so, u must have the form f(r)e where 


mir) 1S approximately constant for large r, and (7) reduces to 
| 1 2 
fern — Oi yee ( tel ye) en © 
v f 
Prrgmror farge r, if f is nearly constant f£ << kf and this equation 
Can be approximately integrated, giving 
; Ps 
2ik mm f = \Cav(rys pas heey A 5 Jdr 
and the right hand side tends to a constant for large r only if 
Meeyegecs to zero faster than l/r (ref. 22, p. 23). 


In the region where V(r) can be neglected equation (7) reduces 


to the spherical Bessel equation whose solutions are 


uy (rT) 





= Ay jy (kr )+ By ny (kr ) (8) 


ee. 





where 3, (kx) and n, (kr) are the spherical Bessel and Neumann 


functions, respectively (23] whose asymptotic forms are 


Dg ) coe Sin(kr-£1/2)/kr 


ny(kr)==3 cos(kr-£1/2)/kr . 


Spnee ny (kr) 1s not well behaved near the origin, if V(r) = 0 
everywhere then B, mast be chesen to be Zero in order for vine 
solution Uy to be valid. Thus it can be inferred that the value 
of By compared with Ay will be a measure of the strength of the 


potential and hence of the scattering. To make this explicit (8) 


can be rewritten, using the asymptotic forms above 


A B 
L | Ln £ Qn | 
uy, (x =a Sin(kr- B= A, cos(kr- 5) 
and letting B./A,= - tan6, 
kn 
DYG) oo Cp nt 64) eo) 


so that equation (6) becomes for large r, a plane wave shifted in 
phase from that of an unperturbed (V=0O) plane wave. Cy can be 
taken equal to one and the normalization absorbed by the expansion 
coefficients in equation (6). 

Tiemttrst term on the right hand side of (1) can be expanded 


as follows: 


ikrcosO _ - lies 
e = Pee i jy (kr )P, (cos®) GO) 


See 


2 ; 
The Neumann function goes as 





a 
=, bel” he] 





and letting Jy, assume its asymptotic form, (6) and (1) can be 
equated and the result solved for £(98), giving, on choosing the ay 


so aS tO assure an OutgOing wave, 


16 


£(9) = (2£+1)e “sinb )P, (cose) (11) 


which is the usual form given for the scattering amplitude in terms 
of the phase shift. 

Deferring the calculation of 7) to the next sections, the sub- 
ject of how the above formalism must be modified in the presence 
of a long range coulomb potential will now be considered. Recall 
that the classical solution for the motion of a particle in an 
inverse square central force field is a hyperbola. According to 
an argument due to Gordon (ref. 22, p. 54), if one considers the 
family of classical hyperbolic trajectories with one asymptote 
originating at the left limit of the z axis, the surface perpen- 


dicular to these hyperbolae goes as 
1 


Z+ am k(r-z) = const 


for large r. Thus it is as if the incident wave were initially 
distorted by the presence of the scattering center at infinity and 
its expected form would then be 


t 
ikl z+ a Ue aera 


a 


A similar argument with regard to the scattered wave leads to the 


form for the asymptotic wave function in the presence of a coulomb 


- 


field 


n aye 1a7 2kr ; (12) 


8 


¥(r) oiket lan k(r-z) 


Yr @ 


iS 





The coefficient a will be defined below. In solving equation (7) 
in the presence of the Coulomb potential one typically converts to 
parabolic cylindrical coordinates whereupon the Schr oedinger 
equation becomes a version of the hypergeometric equation whose 
Solitdenscesare the hypergeometric functions (ref. 22, p. 5/7 ff). 
Using the asympototic forms of the solution functions as before, 


equation (9) becomes, in the case of the Coulomb potential 
xr) —> C,Sin(kr zi ree IG + + 8 if 
where 6, is as defined above and 
Hpe= arg (1 (241 eeia)) (14) 


is the Rehee Sa@ete at Jnismity cue to the Coulomb potential. 
Equation (11) must now be modified to account for this additional 
phase shift (ref. 22, p. 65 ff). 

@re coerficient a in (12), (13) and (14) is proportional to 
the charge of the scattering center. However, Since the problem 
at hand is one in which the pure Coulomb field exists only at long 
range, this charge is the net charge of the ion, in this case (Z-1) 
where Z is the atomic number of the scattering center, and in the 
system of units in use a can be written 

a =IGZ=P)7k (laa) 
and now 5, can be interpreted as the phase shift due to the de- 
Peeture iIrem pure coulomb scattering at close range. It is this 
quantity which is of physical interest. 

To summarize, it has been shown how the asymptotic form of the 
wave function (equations 1 or 12) can be calculated by solving the 


Schroedinger equation for a series of linearly independent functions 


16 





(equations 9 or 13), in other words, that the problem can be reduced 
to solving a series of "partial wave" equations for the various 
values of the angular momentum quantum number £. A simple kine- 
matic argument indicates that as the energy of the incident beam 

is reduced the higher order terms in the expansion (6) become less 
important, enabling a fairly good representation of f(9) to be 


constructed uSing only the first few terms. 


B. THE KOHN VARIATIONAL METHOD 

Before taking up the Harris method it will be worth while to 
briefly examine the Kohn method since it is the most widely used 
variational method in low energy atomic scattering and its dif- 
ficulties were the motivation for Harris' development. These dif- 
ficulties have been examined elsewhere ho29 | and will only be 
mentioned here. 

The Kohn method provides an estimate of the phase shift which 
iow sCeond GMicr accurate’ in the error, i.e., the error term is at 
least second order small. In developing the Kohn method one assumes 
a trial function of the form 


© = sin kr + t cos kr + ¥ (Cicn 


where t 1S an estimate of the tangent of the phase shift and Y is 


normally a linear combination of some basis set X such that 


Y= a.x. — 20 (17) 
1 


So that 1S now a function of Nt+l parameters--the N values of Ass 
and t. Let YW be the true wave function whose asymptotic behavior is 


WV eegeea) Kt tetandscosakr 


W(0)= (0) = 0 
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and define the error wave function 
e=o -W, (18) 
Now form the functional 
F = \ ~(H-E)p dr (19) 
which, uSing the fact that (H-E) #0, can be transformed into 
F = \ P(H-E)€ dr . 
Integrating twice by parts yields 
1 r — 
Fo=- Be - €VD)+ ( E(H-E)@ dr 


pa@esubstituting (16) for € in the first two terms and for 9 under 


the integral gives 


=> k _ 

eS \ M(H-E)® dr = 5(t-tand)+ \ E(H-E)€ dr 
or 

k k 7m 

5 tan6 = 5 t= ho \ S(ier )e- dr (20) 
which identity 1S due to Kato [24] and exhibits the nature of the 
error term. Now the Kohn prescription takes as its functional 


k _ 
ES < t - \ O(H-E)Y dr (2m) 


and seeks its stationary value subject to the conditions 


OT ont ee 
a ys Os 1=l, SN 


She e-alue Of TI thus found is the estimate of tand to second order. 


Co. THE HARRIS EXPANSION METHOD 


1. Formal Solution 


It can be shown that the development of expressions 


for the phase shift lose no generality when the system is restricted 


138 





‘fe one of s-waves (£=0) scattered from a central potential. There- 
fore this restriction will be imposed on the present derivation of 
Harris' method [10] and the Walagdity of its extension to atomic 
systems: will be assumed. 

With the above considerations in mind, assume a system 
describable by a potential V(r) and a Hamiltonian H. Construct a 
trial basis space using a set of N linearly independent functions 
xe. In this restricted space the Hamiltonian can be represented 


by a matrix whose elements are 


es <x bulx > (2) 


while the inner products of the basis vectors are 





ie Se. 
Le X 5 


a 2 
- Ss, (23) 


Then the eigenvalue equation 


N! 
© 


(H - AL)C 


mete the elements of H and L are as given in (22) and (23), can 


be solved for the N eigenvalues A, and corresponding eigenvectors 


p 


Cy2(P=1,..5N). From the eigenvectors a new basis set is con- 
Structed (which spans the same space as the Bea) 
. N 
lo,> = Bip Xs? , (25) 


Now if the true wave function is represented by 
W = s+tc + $ 
where V(0)= 0 


and : 
ee) pas 10 Piet tano 1G¢Os kre — SG ttc 


ee 





(where now the identification of the terms S, t and C is Clear) 
Andent in addition | > may be reasonably well represented at some 


energy E by |§7, a linear combination of the lon? 


N 
Is> = 5 blo > 


then the approximate wave function may nearly satisfy the 
Schroedinger equation. If this is to be the case, then in the 


space spanned by the one the vector 


(H-E)|S+ tc+ & (26) 


must be the null vector, which 18, of course, merely the statement 
of Schroedinger's equation in the restricted ee of the trial 
fiimetilons. This condition may be satisfied if it is required that 
(26) have no component in the space spanned by the Ke that is, 
iE hact 


Soni (H-b)|Seate+ = = O, Ova. 27) 


QO 


Since the lp > are themselves linear combinations of the original 
basis set [ees the condition (27) is equivalent to requiring that 
(26) have no component in the space spanned by the [ex Bearing 
Piem@ind that the solution sought involves finding t, an estimate of 
tan6, note that if <o, | H-E]§? and <p, |H-E|S+ tC? are separately 
equal to zero then (27) is satisfied and t can be immediately found. 
Note further that if E = do the eigenvalue found in (24) corre- 


sponding to the vector 0 that then 


<@,| (H-B)|§> = 0 (28) 


20 





for 


N N N 
- >i ST : ae po yes 
epi (H-E)IS? = Eby ee, XT UALIIX > Cap 


and the third sum in this expression is merely (24) in component 


notation, whence 


N 
: ~ a tee. = (OR ati Cerra 
Mae (HAA )IX Cio j N 


irrespective of the value of Ww. Hence (28) is satisfied and it 


follows immediately from (27) that 


<o | (H-E smc ce 
fano = t = = =—— (29) 


< Hee ike 
oda p)| 


and the problem is formally solved. Finally as the size of the 
basis set ee 1S increased (i.e., approaches a complete set of 
Quadratically integrable functions) [E> Should more closely 
approach |? >and t should converge to the correct result. 

Kolker [25] has shown that the Harris method can be for- 
mulated aS a variational principle and that under certain con- 
ditions it can be combined with the Kohn method to give a minimum 
principle. The advantage of the Harris method 1S principally con- 
putational, in —_—— integrals of the form 26\(Her)ie= need be 
evaluated. However the contribution of these integrals is needed 
to provide the second order correction to the estimate of tan6 in 
variational methods dependent upon the Kato identity (20), so it 
can be seen that the Harris method is necessarily less accurate 
than the Kohn method, as was found by Nesbet [8]. Ilteremanns truc, 


however, that the correction term can be made arbitrarily small by 


all 





wmemeasuna the Size of the trial basis set, for as this set 
approaches completeness, the error must approach zero. 
Qa epplieation to Atomic Scattering 
In applying the Harris formalism to the problem of atomic 


scattering the Hamiltonian is written as follows 





lie 2 are 
i Za ee (30) 


where the subscripts refer to the individual electrons, Z is the 


=? ob 


atomic number of the scattering atom and Tio 7 15, 7 Ts 1s the 
inter-electron distance. The zero of total energy in this system 
occurs with all three particles at rest and separated from one 


another by an infinite distance. The trial basis set is taken as 


the Hylleraas-type functions [20] 
Of 
cs ataset, eo Hered) 
x. = Ge + x 2)y™ e 2 mal ne (31) 
2 
where n, £ and m are integers, chosen subject to the constraint 


n+ +mSmM, (M=1., 2.0ee® ) (32) 


and @ 1S a variable parameter. 

This form of the wave function allows for electron exchange, 
l.e., the interchange of roles between bound and free during the 
scattering process, and enables the wave function to be properly 
Symmetrized when the electron spins are anti-parallel (spin = 0, 
or "singlet" state) or anti-symmetrized when the spins are parallel 
(2pem - 1, or “triplet” state) thus satisfying the requirements of 
the Pauli principle F2GN27) This is the only way the spin enters 
the problem and so once the symmetry of the wave function is prop- 


erly chosen no further concern need be taken with the spin 


Ze 





coordinates of the electrons. Inclusion of the Lio term in the 
potential and in the wave function takes account of the induced 
polarization of the atom arising from the repulsion felt by the 
bound electron as the incident electron approaches. Finally, the 
exponential term forces the functions to zero as tr, or ry increases 
so that all boundary conditions required for these functions in the 
last section are met. As will be seen below the parameter @ is 
varied by inspection so as to optimize the results obtained. 

The number of functions in the basis set is found by taking 
all possible combinations of n, £, m subject to the constraint (32) 
for a given M, and eliminating all terms which duplicate the 


function or require it to be identically zero irrespective of the 


emeice of Sign. The results for 1SMS8 are shown in Table I. 


— 
ce iy 


No. termS in basis set 


Singlet Trapillet 
1 3 i 
Z 7 3 
3 Les: 7 
i, oD 13 
. 34 22 
6 50 34 
7 70 50 
8 95 70 
Table I. Basis Set size for various M. 


The choice of coordinate system is of importance in sim- 
plifying the computational details. For the present computation 
the optimum choice seemed to be the one in which the principle 
coordinates are the radius vectors of the two electrons and the 
angle formed by the radius vectors. The six coordinates governing 


motion of the center of mass and orientation of the system in space 
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are cyclic, and the Hamiltonian becomes [28] 








Meo eo 2 ffi 0 2. | a 
i = 2 a ind_. Oe Son 
ar ror, dr- = as r re aD eZ le 06,5 
1 2 ] Z 
_ 2a — (33) 
a ow 12 
her be - (x? + ae = 2) sax cos@g )? by the law of cosines and 
ee 2 1 2 ew. 2 y , 


the volume element for intearals in this system is 


on » 


ie eae 2 ; 
= sfis 
dr, dr, § rj dr jxodrsing, ,de,, (34) 


The advantage of this system is that all of the integrals arising 
from forming the matrices H and L (equations 22-24) can be done 
exactly, while those arising from <p | (H-E)|s+tc> can be done 
exactly in the case of hydrogen and reduce to single integrals 
which then can be done numerically in the case of helium ions. 
Because the electrons are identical and thus formally 
indistinguishable it is not necessary to integrate over the entire 


first quadrant of the r I. plane. It is sufficient to cover only 


1? 
the region t, =1, {29}. To see why this 1S true one must examine 
the integrals generated in forming Hy and Ley using (31) and (33). 
Carrying out the indicated operations gives rise to 52 terms for 


each ae and four terms for each ae Each term consists of an 


integral of the form 


iS Wp eee 
aa S 12 ee V2 Nee ae : 
A(v,A,u) = \ar ar, e re me ee (33) 
Or 
=JO (Guat) 
4s ee 2 V-2 4-2 u-2 
B(V,A ,u) = \ar a, = Yr, 5 Ti5 cos®,,, ; (36) 
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The complete expressions for ee and Li are given in Appendix A. 
Examination of the terms comprising ie and L. Bodicates 
that every term involving A(vV,A,0) and B(v,A,0) always has a 
@eertservent Of zero, Since the angular integration of B(v,A ,2) 
vanishes, B(v,A,2) = 0, so using the recursion relations proved 
in Appendix B it is sufficient to explicity evaluate only the 
A(v,A,2) in order to find all the terms required to generate the 


heeeeand |... .. 
oui leg 


Setting & = 2 in (35) and performing the angular integration 
over the range 0 toT, gives, in terms of the integration region 


outlined in Figure 1 below and the volume element (34) 


A(VSA52) = 161? \ rv r, e Gr) = Ga) 


O 


Inspection of ie Shows that for each term of the form (37) there 


ms amethner of the form 
: og) ©) eer 
A(A,V,2) = 16.0? ( dri, e ,' anlar (38) 
O 


which 1S the mathematical equivalent of integrating (37) in the 
region r ., = ¥, as can be seen immediately by making the variable 


change r, 2 r, in (38). Thus as a result of the exchange of 
particles possible because the electrons are identical, the in- 
tegrals (37) are sufficient to cover the entire quadrant. This 
argument 1S extended below to the integrals arising a 


=p | (H=E)|S+tc>.. 


Equation (37) may be integrated directly, giving 





» 
A(vV,A,2) = =~ I] vy! - =F Ss 7) 
po 2 =O eee 
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Figure 1. Integration Region for A(V,A,2). 


The factor 16 1” appears in every integral; it can be dropped from 
the calculation since it will merely give rise to the same constant 
feeeor in both Hh and _ and will not affect the solution of the 
problem. Because of the way the recursion relations are structured 
the factor 1/@ always occurs in the form wee and thus can be 
factored out of the integration and incorporated into the coeffi- 
cient of each term in the expression for a and Lae Simee this 
is the only way the parameter @ enters these integrations, evalu- 
ating all the integrals need only be done once at the beginning of 
the calculation for a given basis and the values thus found can be 
used over again as @ is varied. 

Unfortunately a similar simplification is not available with 


the integrals from <_ |(H-E)|S+tC>. However as will be seen only 


fe) 
a relatively few integrals must be done numerically for each value 


of @ and k and so the time required to complete the numerical 
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integration (which dominated the time required to evaluate the 
helium ion phase shifts in all but the largest basis sets) was kept 
within reasonable bounds. Since the present calculation 1S con- 
cerned only with the s partial wave (£=0) further discussion will . 
assume £ has been set equal to zero in all equations of the form 
(9), (13) or (14). Furthermore, since the proper form for hydrogen 
results if (13) is taken as the asymptotic form of the wave function 
and Z is set equal to one in (15), the remainder of the development 
will be concerned exclusively with establishing the Harris method 
in the presence of a long range Coulomb field. 

Bearing these comments in mind, note that as r goes to 
zero (13) does not behave as well as one would like. Because of 
the presence of the logarithmic term in the argument of the Sine, 
aS Y approaches zero the interval between zeroes of the function 
approaches zero and u(r) does not approach a definite limit. To 
force u(r) to zero a factor is normally chosen which goes to zero 
fast enough as r goes to zero to override the effects of any 
Singularity which may occur, and goes to unity as r gets large in 
order not to affect the asymptotic behavior of the function. For 
the present calculation the factor was chosen to be (eee /2))- 
This choice was made to insure that u(r) would go to zero fast 
enough that its oscillatory behavior near the origin would not affect 
the accuracy of the numerical integrations involving the asymptotic 
part of the wave functions. 

With these preliminaries aside, the form chosen for the 


asymptotic part of the wave function may now be presented: 


2/ 





2 
2 


- 2X = r 
StuG.= e 1(q-e 2)3| sin( kx, + Z=* On 2kr + =) 


Zr1 )| 
tand cos( kr + ¥ Qn 2kr + n /x ., 





+ 

6 4 
=i Py - =- Y 
2 io I a Tjek 

ae: (1 -c ) | sin( kx + neds 2kr,+* n) 

Zs aves cos( kr a 22 tn 2Knx * n) |/2 (40) 
i kK as: al 

where : 
n= arg {r(1-i 224 )} a2) 


Evaluation of the complex Gamma functions is discussed in Appendix 
D. Normalization factors in (40) have been suppressed because they 
@e not affect the calculation, cancelling when t is evaluated. 
Following the prescription outlined in section 1, the 
estimate of tan6d is evaluated from theratio of the integrals 
(uH-E)|S> and <o@ 


<o gies) fee at the energy Eo found by solving 


p p 
(24). Since the functions (H-E)|X > have been previously evalu- 
ated in solving (22) and (23), considerable simplification is 


achieved by taking advantage of the fact that H is a hermitian 


operator and solving instead the equivalent integrals 


N 
 .* | (11-E,) ae (42) 


and 


N 
a ee | (H-E,)IX, > (43) 


whence t is found from the negative ratio of (42) to (43). Note 


again that 5 measures only the departure from pure Coulomb 
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Peer texans arising in the region close to the target where the 

three particles in the system must be treated explicity. 
Carrying out the indicated algebra one finds that 

<sl(H-E)|\¥> has 52 terms constructed of integrals of the form 


racs e a 3 
eos! eee oo . 


. Z-1 _— = 
sin(kr + i in 2kx,+7)) dr, dr, (44) 
OL Ol Of 3 
ae oo Be 2S 528 (. -271) : 
Bet) i 329 “Te i 
: Z-1 _ <> 
sin(kr,+ : in 2kr, +7) dr, dr, (45) 
OL o OL 3 
ae ie.) = es et 22 ( "2 *2) 
iia lame hee : aS a 
: Z-l1 , 7 
Sin(kr,+ mes 2kr,,+ n)cos®,,dr, dr, (46 ) 


OL OL 3 

= tee ee = = = 
Vee ecee b=? ( 22 2 ( ) 

B,(V,A ju) = \7Y t T)2 e e ! l-e a x 


Sin(kr,+ Z°* tn 2kx + n)cos®,, dr dr, (47) 


Four similar integrals with cos(kr + ....) replacing the sine terms 
in (44) - (47) arise in evaluating < cl (H-E)|x > and are designated 
Ag; le Ba) and By? respectively. Note that.all these integrals 

are of the same form as those in Appendix B Section 1 and hence the 


Same recursion relations may be applied. 
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As with the integrals of equation (34) the evaluation of 
(44) - (47) starts with the case # = 2 and as before, the angular 
integration yields B, (vA 52) = B, (vA 52) = O and introduces a 
feetor 2 in the coefficient of A, (vA 52) and A,(V,A 52). Also as 
with the A(v,A,2), terms of the form A, (vA 52) and A, (vA 52) 
always exist together and thus satisfy the requirements that allow 
the integration to proceed over ry) = I. only, as illustrated in 


Figure l. 


Starting first with A, (v5 52), after angular integration 





(o% r e4 (eG 3 
ace -(Z+ =)x 1 — ( Scar ) 
2 ? 
A (924 52)= 167 \ dr re s a dr oe 2 2 l-e ae Os 
B ceed Ze 
O O 
Sin(kr,+ 2-5 on 2kr + N)- 


Now referring to Figure 2, note that the order of integration can 


be reversed by choosing the elements of area dr, dr, as indicated 


while the same region is covered. With this change Aj, (vA 52) 


becomes 


7 O Ox 3 
2(° ra * ol — 3 ey 
1 rene 


— a 
Sin(kr + on On 2kr + n)\ ced cl ea 
r 


2 





and now the integral over r, can be done giving 


i 











S 
ek. cya te gore, a? Ss BD) 
Av; »2)= 71 ye Ge dr 5t., l-e 
7 ee) 
é Z-1 i) 
x sin(kr., ce vasa dd 2kr.+ 19) 
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Finally, making the variable change x = (Z + d)r,, 





Vv 
320%! BZ Zo iL 
Bye 22) = VtA 4+] |e (322) (v-i)! i aoe co 
(2Z+0) (Z+0) 1=0 . 
where 
60 4 - Basa). 
a ee = (ZH e& Z-1 yp, 2kx ) ; 
1,(Y) \.° So heres sin a40r te 7” ie = iG (49) 
* . 
when Z =2 (He ), equation (49) must be evaluated numerically. How- 
Svetrmichn Z = 1,7 = 0 and the integral reduces to a standard form 


(30). Numerical evaluation of (49) is discussed in Appendix D. 


bg 


2 
© co 
1 foraf ary 
O r 
rae ~ 
Z, 
4% 
2 
a 1 
| fy f dr, 
O O 
I 
Figure 2. Interchanging the ry and Ir, Integration. 


Proceeding now to the evaluation of A,(V,A 52), note that 
the r, integration may be done directly and that evaluating it at 


the upper limit (r,) Gives the same function of r, as did the 


evaluation of the ry integration at the lower limit (r.) in the 


evaluation of Aj (vA 42) but with v and X interchanged: 


gill 











L Bees, ! (52 )\ (Goage 1) 2 im. ull 
Aj(V5A52)= 274+ a 3740 dr l-e e rj sin(kr + — 2kx 5 +7) 


y a a \ ( “3") ee even Z-1 

= = : - Sen =n 7 } 
(32 RE Po ote eC ry sin(kr,+ 7 ekr, n ) 
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Making the variable change x = ar, /2 ijeathieoeticstetermiand 
x = (Z+Q)x, in the remaining terms, and recognizing that the 


summation is in fact A, (AV52) gives 





ret 
Ao(V5%,2)= _32n" = ——+ T_(v)- A, (A,v,2) (50) 
(2Z2+2)a& (27405) 
where 
a \= A aa 14 i “xX 3sin( 2 + Z-1 wn 4kx os oF (51) 
aah 4 x ( e ) > 1 OY k ; Of 1) se 


Now since A, (v,A,2) and Aj(A 532) always occur together it is con- 


venient to define 


A_(A ,V32) = A, (v.A,2) + A(A,Vv,2) (652) 
and similarly 

A (A,v,u) = A, (vA jt) + Aj(A,V4H) (2329 

B.(A,V,4) = B) (vA 51) + Bo(A, Vu) (535) 

JANE) A3(V,A,u) + Ay (A svat) (53¢) 

B(As vob) = Bz(vA 5) ey By (A 5v 51) (530) 


where (53a) and (53b) arise in evaluating < §|(H-E)|x> and (53c) and 
(53d) arise in evaluating <c] (H-B)|x > and the number of terms in 
each is now reduced to 26. Unfortunately, explicit expressions for 


A, (v,A 52) and A,(V 5A 52) are required in order to evaluate A.(vA,1), 


Be 





B.(v,A,1); A (vA 51) and B.(vsA 51); respectively so the reduction 
in labor is not so great as it might appear. Note that the factor 


32n*/(2z+a) is common to every term, so it may be dropped from the 


calculation. 
3. Summary of the Harris Method 


To summarize the Harris prescription as it applies to 
atomic scattering, the computational scheme is as follows: 

a. For the given basis set Size evaluate the required 
Pelose-im” ingegrals from * 


r i 
a , di 2 (Atv-i)! 
A(V,A 2) or a {v: = es ae iB =) } (54) 


b. Evaluate the required A(v,A,uU) and B(v,A,u) using the 
recursion relations of Appendix B. 

©. For a given @, evaluate ig and at using the expressions 
: ee Neca! 
in Appendix A. Note that each A and B must be divided by @ 
mo give the correct result. 

d. Solve the eigenvalue problem (24) for the eigenvalues 
ES and corresponding eigenvectors ah 


e. Find an appropriate incident electron momentum from 


the expression 


k = 2E +2 (55) 


Cees 
where -Z gives the ground state energy of the scattering center 
mi rydbergs (1 ryd = 13.6 eV). 

f£. For this value of k and the previously assigned value 


of @, evaluate (numerically if necessary) I Dijana vee for the 


1’ “s 3 


required values of the index Y using (49) and (51) for I, and I, and 


il 





> ————~ 3 
7 y( 2(Z+a) ) ( Kee 7-1, 2kx ) 
T,(¥)= = pra a cos\s ae! oa oP gay iebs (56) 
and 
™ y 3 k t,t 
= a. a ee Sa fp ) 
ey) ‘ ae e .aeos ( - + z Yn 7 + 9} dx (ow) 
g. Using these values, evaluate the required terms of os 


A,» A, and Ay for tw = 2 using (48) and (50). 
h. Using the recursion relations in Appendix B evaluate 
the remaining required terms of ee Bos ave and Bo: 
1. Evaluate the < Sy (E18) aes and < c| (H-E)|x, > using 
1 
the expression in Appendix A and form : c<S| (H-E) |X, 7 and 


1=1] 


N 
7 ee (4eR) be 5. 


j. Finally the negative quotient of the two sums above is 
ey the desired result, and repetition of the above step for 
various @ normally permits the entire desired range of k_ to be 


covered. 
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iii hes ul Scans DISCUSSION 


The calculations reported in this section were performed in 
double-precision (16 figures) arithmetic on the IBM 360/67 computer 
of the Naval Postgraduate School Computer Facility. The programming 
language was Fortran IV and the H-level compiler was used. The sub- 
mootanes used to factor the matrix L and invert U and to solve the 
eigenvalue problem, and the 32-point GaussS-Legendre quadrature 
Subroutine, aS well as Several auxiliary subroutines were furnished 
by the source library of the Computer Facility. 

The principal program was, of course, written to evaluate the 
phase shifts. This program and its required subroutines are re- 
produced following the appendices. In addition it was found con- 
venient to have available modifications of the main program which 
solved only the eigenvalue problem, giving in one case the values 
of the incident electron momenta in the elastic scattering range 
and in the other the calculated energy levels of the bound state 
system. Since the scheme adopted for checking the eigenvalues and 
eigenvectors was rather time consuming and involved a rather large 
quantity of printed output, it was done on only selected matrices 
and then only in the supplementary programs. The modification of 
the phase shift program to perform these supplementary functions 15 
obvious and the programs so modified are not reproduced here. 
Several features were built into the programs to permit the maximum 
flexibility with the minimum of internal changes and as it finally 


evolved only the storage capacity within the program had to be 


315, 





changed, depending on the size of the basis set to be used. All 

other options were controlled by data inputs from external devices. 
The results for e-H scattering will be presented first, followed 

in part B by the e-He” results. Comparison with the results of 

other calculations will be made where appropriate. The results 

will be presented in the following order: (1) phase shifts for 


Singlet and triplet scattering; (11) location of resonance levels; 


and (111) bound state energy levels from the eigenvalue problem 


(see Appendix C). 


A. ELECTRON-HYDROGEN SCATTERING 
1. Phase Shifts 

Singlet and triplet phase shifts were calculated for basis 
sets from N = 22 elements (M = 4, singlet; M = 5, triplet) to 
N = 70 elements (M = 7, Singlet; M = 8, triplet). Although the 
program had built into it the capability of calculating phase 
shifts for the 95-element basis set (M = 8, singlet) the eigen- 
value problem began to show signs of instability at this size and 
the computer time required to find eigenvalues and eigenvectors 
became unacceptably long (8-10 minutes per matrix), permitting no 
more than testing the program at this level. In addition the 
storage area required for such programs approached the limit avail- 
able on the computer. 

Preliminary to calculating the phase shifts, sauereion ( 24) 
was Solved for a wide range of values of the non-linear parameter & 
and the resulting values of incident electron wave number were 


plotted as a function of @. As expected, the eigenvalues fell into 
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Sy discernable trajectories whose evolution as N was increased 
was eaSlily traceable, at least in those regions where the eigen- 
values were sparse. Figure 3 shows such a plot for the singlet 
states of hydrogen with N = 70. 

This plot clearly illustrates a major disadvantage of the 
Harris method, which has been alluded to by others [13,25] in the 
context of the inconvenience of being unable to pick the scattering 
energy without the expenditure of considerable labor. The dis- 
advantage seems to be of a more fundamental nature, however. While 
coverage of the entire energy range 1S certainly possible, at least 
in most cases (it was not possible with the basis sets used to reach 
zero energy in the e-H triplet case), sufficient points at any given 
value of k are not always available to permit the sort of investi- 
gation of convergence rate that characterized the analyses of 
Schwartz [aa and Armstead eanle Furthermore in those regions 
where few eigenvalues exist there 1S no guarantee that they will 
exist for an optimum value of @. Schwartz has indicated that the 
useable range of @ was between about 0.8 and 4.0 in his calculation, 
and this generally proved true in the present case also. Reference 
to Figure 3 shows that even in the largest basis set used the value 
of @ for the only trajectory passing k = 0.1 is around 5.0, just 
outside this range. In thiS region, increasing N is not likely to 
produce additional eigenvalues at which to solve the phase shift 
problem, for as 1s shown in Appendix C, the eigenvalues result from 
a Rayleigh-Ritz variational calculation of the energy levels on the 
two-electron bound system (see part 3 of this Section) and hence 


the eigenvalue trajectories shown in Figure 3 represent approximations 
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Figure 3. Eigenvalue Trajectories for Singlet e-H System N 
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to the various energy levels of H for a given @. Although the 
scale is not the energy scale one uSually associates with such 
plots the relationship between the two is 


wae sol (58) 


in the case of hydrogen, where E 1S given in atomic units (27.2 eV). 
In such calculations, as N increases the lowest lying trajectories 
become quite good approximations to the energies of the lowest lying 
eigenstates of the bound system [32] after which further increase 
in N will‘ give rise to no additional trajectories in that region. 
In fact the situation is likely to get worse from the point of view 
of the scattering calculation, Since as N is increased the minima 
in the trajectories become broader, approaching horizontal straight 
lines as completeness 1S approached. Consider an incident electron 
wave number near but above that corresponding to the highest well 
defined energy level of the bound system for a given N. As N is 
further increased the values of @ for which such a wave number is 
an eigenstate of the system will be pushed closer to, and perhaps 
beyond, the limits of therange in @ for which reliable results can 
be expected. Figure 4 illustrates this behavior for a typical 
elgenvalue trajectory. AS N 1s increased the curve becomes broader 
and its minimum approaches a limit. For N 7 70 the minimum probably 
will not decrease appreciably but the curve will become broader. 

One might conclude that in an energy region which is rich 
in elgenvalues the situation ought to improve. That such is not 
always the case can be seen by examining Figures 5 and 6, showing 
singlet phase shifts for k = 0.4 and 0.8 plotted as a function of @. 


It must be noted that the lines connecting the points arising from 
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Figure 4. Typical Evolution of One Eigenvalue Trajectory as N 
1s Increased. 
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the calculations for a given basis set are just that and no more. 
They have no physical or mathematical meaning, for the points 
plotted represent a complete collection of the possible values 

to be found at a given energy (possibly neglecting one or two 

at small @). 

In Figure 5, while the "curves'’ of all basis sets have 
clear cut minima, there is no clearly defined trend as the basis 
set size Pies. as Schwartz and Armstead observed in their Kohn 
calculations. This may be due to numerical errors which accumulate 
more rapidly as the basis set is increased in size, however as will 
be seen below the ambiguity 1S considerably reduced in the triplet 
case which leads one to conclude that this error is probably small. 
Note however that the range of @ over which the phase shift 1s 
nearly constant 1s greatest for the case N = 70 and hence the 
minimum of that curve has teen selected as the most probable value 
in this calculation. The rather large error indications in Table 
II reflect this situation. 

Referring now to Figure 6, note that while the minimum of 
the N = 70 curve is at 6 = 0.761, there is a plateau forming 
around 6 = 0.886. For this value of k increasing N may well bring 
improved results since the minima of more trajectories may be ex- 
pected to move into the region, however the ambiguity will remain 
if the plateau becomes broader and stable in value. A similar 
Sfructume iS evident in the Galculations for k = 0.7 and 0.6 as well. 
It 1S interesting to note that these plateaus correspond to the 
values for the phase shift given by Schwartz and the difference 
between these plateaus and the minima of the curves increases with 


increasing k until at k = 0.866 the difference is in excess of 10%. 
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Figure 5. Singlet Phase Shift for e-H as a Function of @ 
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Although they do not discuss it, close examination of the results 
of the Harris calculation of e-H singlet phase shifts by Oberoi 
and Callaway [13] using the close-coupling approximation seems to 
indicate that their results may be as much as 4-5% less than 
Schwartz's at k = 0.8, although their graphical presentation of 
results makes such a comparison difficult. 

Appeal to experiment to resolve the ambiguity is fruitless 
at this time for the difference in total scattering cross-section 
at k = 0.8 for the two values of 6, 0.886 and 0.761 is only about 
3%, far less than the accuracy of any experiments conducted to 
date [33]. The experiments report results about 15% less than 
those of Schwartz and Armstead at the upper energy limit and it is 
worth noting that the aera t result 1S also lower than that of 
Schwartz's. However because of the large experimental errors re- 
ported (2 20%) the significance of this fact, if any, cannot be 
estimated. 

Based on the results available it seems more wishful than 
logical to ascribe more significance to the plateaus in the N = 70 
curves than to the minima and in the absence of the results of 
Schwartz the lower value would most likely have been chosen as the 
most probable value. Therefore the minima will be chosen for 
consistency. 

The extrema of the curves for phase shift vs. @ generally 
are minima except at the lowest energies calculated and they seem 
to trend downward, so in most cases the value of the phase shift 
for N = 70 can probably be taken as an upper limit except for 


k S 0.2 where it seems to be a lower limit. Where an unambiguous 
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convergence could be established it Seemed to be about as the 


sequence 1/N. 


The scattering length is defined to be 


tand 
Oo 


ani es aa ‘ (59) 
k70 


In the present calculation 1t seemed natural to evaluate the 
scattering length for e-H scattering by simply evaluating (tané)/k 
for a succession of small values of k and extrapolating the results 
to k = 0. This was done for the case N = 70 and the smallest value 


of k used was .000412. The extrapolation to k = O gave 
a. = ~6.02181(1) 


where the figure in parentheses is the uncertainty in the last 
digit. This value differs by about 1% from the value ale -5.965 
reported by Schwartz (12). Note that the uncertainty reported 
abpewe 1S the uncertainty arising from the extrapolation to zero 
and is not necessarily an indication of absolute accuracy. The 
extrapolation used in evaluating ao 1S shown in Figure 7. 

It should be noted that the difficulties with the 
Harris method discussed above should not be peculiar to this choice 
of trial wave function, for since the Hamiltonian (30) 1s the exact 
Hamiltonian for any two-electron system the eigenvalue calculation 
Should lead to the same effects described above for any choice of 
trial wave function as long as it meets the uSual boundary condi- 
tions and is based on a quadratically integrable basis. Perhaps 
by uSing some approximate Hamiltonian the eigenvalue situation 
would be improved but the loss of accuracy from this approximation 


would probably offset any gains from improved eigenvalue behavior. 
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Evaluation of e-H Singlet Scattering Length. 


N=70 
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L6 








The results for the triplet phase shifts are much as has 
been described above for the singlet case except that the curves 
seem to trend upward toward maxima for values of k larger than was 
med in the Singlet case. Since in the triplet case the electron 
Spins are parallel the effective repulsion between the electrons 
due to the Pauli principle tends to improve the convergence. Thus, 
While the qualitative picture is Similar to that of the singlet 
scattering, the ambiguities are much less severe. Figure 8 shows 
the eigenvalue trajectories for N = 70. Note that there 1s no tra- 
jectory that reaches k = 0, the lowest one just barely passing 
m= 0.1. Henec it Was not possible to evaluate the triplet scat- 
tering length as was done for the singlet case. Figures 9 and 10 
show plots of phase shift vs. @ for k = 0.4 and 0.8. Finally, the 
curve of 6 vs. k for both singlet and triplet scattering 1s shown 
in Figure 11. The resonances shown on the plot are discussed in 
the next section. In each case one resonance which has been re- 
ported elsewhere as being just below the excitation threshold was 
observed in this calculation to be just above and is so indicated 
on the plot by the vertical dashed lines. It is believed that in- 
creasing the basis set size would bring these values below the 
threshold. The singlet and triplet phase shifts are tabulated and 
compared with those of Schwartz l12} in Tableers. 

2. Scattering Resonances 

The present calculation offers a feature not found in 
Schwartz's calculation--the localization of real scattering re- 
sonances. There 1S probably no theoretical reason why the Kohn 


method should not show the resonances, however since they are very 
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narrow (18,34,35] iris quite likely that they simply were not 
resolveade’ It might also be difficult to distinguish them from the 


spurious Ones Characteristiec™or such calculations. 


Incident 
Electron Scaipmeringe Phase Sait in Radians 
Wave Number Singlet Triplet 
Present Present 
Work Schwartz[12] Work Schwartz(12] 
O.1 aeib6O Sie 2553 2.9386 (5) 2.9388 
OR2 2.067(4) 2.0673 2078S fey) 2 gs. 
0.3 1.696(6) 1.6964 225013 7) 2.4496 
0.4 1.414(6) 1.4146 292930(.5) 2.2938 
Owe 5 1.194(7) 1202 2.1093(3) 2.1046 
0.6 1.020(4) LAO 193215. 1.9329 
0.7 .878(3) . 930 Leo ie) Me oe. 
0.8 .761(4) 886 1.6 389(5) 17653 
0.9 .693(7) --- ke 5622.0) EPs 


a : : : : 
The figures in parentheses represent the uncertainty in the last 
digit reported. 


Table II. Singlet and Triplet Phase Shifts for e-H 
Elastic Scattering. 

The existence of scattering resonances has long been 
associated with the auto-ionizing states of the compound system { 36 ] 
and recent experiments and calculations have confirmed this 
35, 37.38). Auto-l1onizing states (in the case of ae they are more 
properly called auto-detaching) are those in which the total energy 
of the system exceeds the ground state eneroy by more than the 
1onization energy of the two-electron system but in which both 
electrons are in excited states. Such states can decay non- 
radiatively when one electron's excess energy is transferred to 


the other which is then ejected from the atom while the first 


reverts to the one-electron ground state. Such states are normally 


ae 





extremely short lived, typical lifetimes often being several orders 

of magnitude less than the normally allowed radiative transitions [39]. 
Clearly an electron captured momentarily into such a state will be 
considerably delayed in its passage by the atom, suffering a cor- 
respondingly greater phase shift, and thus show a resonant behavior 

in the cross section. Since these states lie above the one electron 
10nization potential and are extremely short-lived they are often 
referred to as quasi-bound (40). 

Figures 3 and 8 show the evidence of such states in the 
fmolet and triplet structure of H . They are the broad shallowly 
curved trajectories lying just below the excitation threshold in 
each case. 

Since these states do lie above the ionization potential 
of H and thus among a continuum of states, one expects from the 
Hylleraas-Undheim theorem [ 32] that such states would not be well 
approximated by an eigenvalue trajectory until after all lower 
lying states had been. The reasons for this mt being so are not 
well understood, however Holgien and Midtal (40) have found that if 
the trial wave functions include an appropriate mix of the one 
electron state functions (i.e., 1S5ns, n 7 1, and also l1sl1s, 2s2s, 
2p2p, ete.) that the eigenvalues corresponding to the autoionizing 
states begin stabilizing long before those of many lower lying 
States. They conjecture that this probably occurs because the 
equal treatment of the two electrons leads to an approximate 
orthogonalization to the lower states. The limits approached by 


these trajectories may not be upper bounds however. 
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That Such conditions occur in the present problem is most 
likely a result of the use of Hylleraas-type trial functions, which 
can be combined to approximate the one-electron orbitals needed to 
satisfy the condition set by Holsien and Midtal above. 

In the present calculation two resonances were resolved 
within the e-H elastic scattering range. One in the singlet states 
was observed at k = 0.83713 which corresponds to a bound state 
energy of E = -.149610 atomic units (1 a.u. = 27.2 eV). A second 
in the triplet states was observed at k = 0.86418 which corresponds 
to E = -.126596 a.u. In addition there were two such states identi- 
fied just above the excitation threshold which probably correspond 
to levels predicted just below the threshold by others [28u, Tn 
all liklihood a larger basis than those used here would give a 
better approach to the results reported by others. Since the means 
of locating these resonances was from a bound state calculation no 
estimate of their width can be given. A characteristic of these 
calculations is that the phase shifts associated with a particular 
eigenvalue trajectory tend to be rather independent. Thus evaluating 
a phase shift at a value of k corresponding to a resonance energy 
but at a value of @ corresponding to a normal (i.e., non-resonant) 
elgenvalue trajectory appears to show little or no evidence of the 
resonance. Although the non-resonant and resonant trajectories 
sometimes do cross, no conscious attempt to evaluate the phase 
shift at such a crossing point was made since in the eigenvalue 
method used (Jacobi variable threshold C41] the eigenvectors found 
for degenerate or nearly degenerate eigenvalues may be grossly in- 


me ctrerve Retr. “1, p. 261) and thus the phasé™shift calculated at 
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this memem woulda be difficult to interpret. The independence of 

mae phase shifts corresponding to different trajectories is probably 
related to the orthogonality properties inherent in the eigenvalue 
problem and discussed above but beyond this it is not well understood. 
Table III summarizes the results for the e-H resonances and compares 
them with other calculations and with such experimental data as is 


available. 


eae eneuay Levels 


As shown in Appendix C and as demonstrated in the last 
Sectron, the first part of the Harris prescription involves what 
amounts to a variational calculation of the bound state energies. 
Reports of such calculations abound in the literature, commencing 
with Hylleraas [20] “ana continuing to the present, and have 
attained a high degree of accuracy. Hence Phe vee obtained 
here are of no more than passing interest except that they arise as 
a natural by-product of the Harris method, and this fact seems not 
to have been mentioned previously in the literature concerning the 
Harris method. Considering that no special pains were taken to 
asSure accuracy (other than the normal and reasonable ones dis- 
cussed in Appendix D) the energy results are Surprisingly accurate, 
and since they are obtained with no additional effort from the 
phase shift calculation it seems worthwhile to discuss them briefly 
here. 

Pekeris [46] has reported a calculation of the ground state 
energies of He and H in which he used matrices of up to order 1078 
and was at great pains to achieve the maximum possible accuracy. 


He reports the ground state of H from this calculation to be 
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pip2775097 a.u. In the present work a value of -.5277477 a.u. was 
obtained with an order 70 matrix. 

As an example of the bound state results obtainable, 
Figures 12 and 13 show the negative energy structure for the 
Singlet (N = 70) and triplet (N = Pa caneuiaioasy Not all the 
detail is included in the region near zero energy because the large 
number of ecigenvalues in that region makes it very difficult to 
sort out the eigenvalue trajectories. Note that the lower eigen- 
values clearly exhibit the broad minima alluded to in section l. 
This characteristic of the trajectories 1S even more pronounced in 


the case of helium, which will now be discussed. 


Se ELEGPRON-HELIUM ION SCATTERING 
i. hase Shifts 

In principle the only difference between Scattering elec- 
trons from helium ions and from hydrogen is the presence of the 
long range Coulomb force in the former. This effect 1S accounted 
for in the form chosen for the asymptotic part of the trial wave 
function and once the more difficult problem of numerically evalu- 
ating the resultant integrals iS solved, one expects the procedure 
to be straightforward and to yield results whose accuracy is de- 
graded only by the relative inaccuracy of numerical integration vs. 
exact integration. Appendix D contains extensive discussion of the 
method of evaluating the numerical integrals and assuring their 
accuracy, and it appears that a high degree of confidence can be 
placed in them. Nevertheless the results of these phase shift 
calculations are not so accurate as the hydrogen ones, if one 


accepts the results of Burke and Taylor [18] as the most accurate 
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Figure 13. Triplet Energy Levels for H .N = 50 
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to date. The close agreement of their results with those of 
others by different methods [35,44] indicates that this confidence 
is probably well placed. Thus 1t would appear that the Coulomb 
force has an effect on the accuracy of the calculation beyond what 
is initially expected. 

As Figures 14 and 17 show, the eigenvalue trajectory situ- 
ation is somewhat improved in the case of helium, owing to the 
existence of a complete set of two electron bound states in the 
region below the first ionization potential of He. Nevertheless 
the singlet phase shifts at low energy appear to be as much as an 
order of magnitude larger than the results reported by Burke and 
Taylor. Because the calculation gives only tan6 the value of 6 
fewuneertain by + n™, and the choice of 6 to be .4 at k = 0.2 is in 
closer agreement with Burke and Taylor. However the slope of the 
6 vs. k curve at k = 0.2 was measured and its value seems to support 
the higher value chosen. At higher energy the convergence appears 
to be to a result about 5-10% higher than that reported by Burke 
and Taylor, although the existence of the scattering resonances 
at higher energies undoubtedly affects the results and makes com- 
parison difficult in this region. The existence of the resonances 
makes more likely the possibility of attempting to evaluate the 
phase shift at a nearly degenerate eigenvalue, with its concomitant 
inaccurate elgenvector which may further confuse the picture. 
ivyprcalmmelots of 6 vs. @ for singlet e-He™ scattering are shown in 
Paoures 15 and 16. 

| A similar but (not unexpectedly) less severe situation 
exists with the triplet scattering. Figures 18 and 19 show typical 


triplet plots of 5 vs. q@. 
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Results of the present calculation are presented and con- 
pared with those of Burke and Taylor in Table IV and plotted in 
Pagure 20. 

Meese attering. Resonances 

It was possible to identify six singlet resonances and four 
triplet resonances. As shown in Table V the agreement with the re- 
sults of Burke and Mevicar [44] and others [37,40,42,43,47] is 
quite good. 

Burke and McVicar give the energies of seven and five re- 
sonances for singlet and triplet states, respectively. As in the 
hydrogen case an additional resonance level could be identified just 
above the excitation threshold by examining the energy level plots 
for each system. Again as in the hydrogen, these levels would 
probably appear just below the excitation threshold if the calcu- 
lation were done with a sufficiently large basis set. 

3. Helium Eneray Levels 

The helium energy Spectrum has been even more extensively 
calculated than has hydrogen (20,46,48,49]. The most accurate such 
Galeulation 1S that of Pekeras [46 ] who has reported a value of 
-2.903724375 a.u. for the ground state energy of He. After applying 
mass polarization and relativistic corrections and correcting for 
Piiemibamb Shatt, he finds a value of 1.98310687 x 10° ona Tore tie 
ground state of He, which compares with the experimental deter- 


5 


mination by Herzberg [50] of 1.983108 x 10° + .05 Se Sa Rey. 


5 


£orerne 2°S state of He, Pekeris obtained -2.17522937822 a.u. or 


3.84566 x 104tem72 while Herzberg measured 3.845473 x 104 + ,05 oe 


— 


In the present calculation -2.90372410 a.u. and -2.175191 a.u. were 
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5caé theming Phase Shift in Radians 


Incident Singlet Triplet 
Blectron Present Bue kesaid — PECSene Burke and 
Wave Number Work Taylor [18] Work Taylor [18] 
0.2 Ce cha aan 

D.4 .98(4) . 1.469(8) 

@.4472 2 7 .9088 
0.6 Bose ) ie LOG to) 

0.6325 ~4092 8873 
6.7746 » 3989 . 8649 
0.8 Seo) .965(6) 

0.8944 3912 .8462 
lees O 446(7) » 3850 .889(4) 8253 
f.1832 3780 .7895 
H.2 Lali 7 } eo SS) 

1.3416 sony? BOIS 
age 407(4) .778(4) 

1.4832 .3910 ec 
B.5611 as 

1.6 348 (4) -730(4) 

-.0553 . 3988 

f.6712 .7006 
6971 .6931 
f .7073 Seo7 2 

ft .7 32 . 326 (3) .705(4) 


wee eee eee eee ee eee ee eee eee 


a : : 
The figures in parentheses represent the uncertainty 1n the last 
Gedit weportced. 


+ 
Singlet and Triplet Phase Shifts for e-He 
Elastic Scattering 


Table SV. 
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Figure 20. Calculated Phase Shifts for e-He Scattering. 
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obtained for the 1° S and 23 


S levels in He, respectively. Figures 
21 and 22 show the calculated energy levels for the singlet and 


triplet S states of He. 
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Figure 22. Triplet Energy Levels of He. N=50 
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Wee CONnCeuUSION 


The principle advantage in calculating phase shifts by the 
Harris method is that it is not required to evaluate matrix 
elements involving only the aSymptotic part of the wave function. 
This results in a considerable computational simplification but 
the price paid is that the error is now of first order and hence 
the method is inherently less accurate than methods dependent on 
mmemiato 1dentity. 

An additional disadvantage is that in certain regions there may 
be a paucity of energy eigenvalues at which to carry out the calcu- 
lation. Since these eigenvalues are approximations to the compound 
system energy levels, in those regions where the approximations are 
good, increasing the number of parameters in the trial wave 
function may simply make the situation worse. Even in regions rich 
in eigenvalues the fact that they are related to the bound states 
may have unpredictable effects on the phase shift calculations. 

The Harris method does have a feature that could prove useful 
in complex systems. Resonances in the scattering cross-section 
can be uniquely identified as arising from certain bound states of 
the compound system, such as the auto-ionizing states in the pre- 
sent instance. Since the scattering at a particular energy can be 
identified with the compound state energy level, when that level is 
well defined by the trial wave function the correlation with the 


behavior of the phase shifts in that region will be clear. 


74 





Whether the anomalous results for the e-He system found here 
are an indication of a defect in the Harris method or due toa 
subtle error in computation or analysis on the part of the author 
memmot absolutely clear. Although the latter cannot be ruled out, 
the success of the method in calculating the e-H scattering using 
the same computer program and even, as pointed out in Appendix D, 
the same integration scheme (on a test basis), Seems strong evi- 
dence in favor of the former conclusions. A third possibility, of 
course, is that the results obtained here are correct and the re- 
sults of Burke and Taylor and others are incorrect. [In the 
absence of experimental evidence to the contrary, the consistency 
of the other calculations would seem to make this conclusion also 
unlikely. 

Thus it 1S the tentative conclusion of this thesis, that the 
Harris method should be applied to electron-ion scattering with 


some degree of caution. 
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APPENDIX A 
Formulae for the Various Matrix Elements 


a H..ands.. 
1j 1) 


The matrix elements < y.July.> is made up of 52 terms, each of 
2 J 


which includes an integral of the form 
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notation a can be written 


Z 
Hs = - (o oe ane n+ 2,4 .+ eee 2,m,+ Bae) 
+ (2 ano “A(n.+ n.t 1,£.+ &. + 
a(n, )-z J on n+ 1,04 Hi 2,m,+ m. 2) 
mae + 1)- . ; : ; ; ; 
L mr ( é )-z J A(n;+ ne 2,4.+ ra 1,m,.+ eee) 
+ } ee : : Meee : 
ame A(n;+ ae Ze a oh 3,m,+ or ) 
ve 1. . £ 
at A(n,+ a 3, ar a 2,m.+ ue ) 
= @nmenet £.t not 1)°Aimat n.t 2,0.¢ £ + 2,m.+ m. ) 
See, J ) 1 J 1 3, 1 s) 


- mete t | 1 \> + ae) 
epee ) A(n; ee i a 2,m.+ Ble 2a) 


76 








+ 


+ 


Uau(4 + Lyez 


Uzu(n + 1)-z] 


<n @ (es 
A(n, : 


+ T).. + 
A(n,; 4 


*B(n,+ ns 
*B(n, + a 
"B(n,+ Nag 
"B(n, + la 
eas ao 
“A(t aru 
ae Da 
eas ree 


; L + 
J cry j 


tien pone ) 
2,4 .+ Se ol “ 


ean 1h) 
2,k.+ a 2,m,+ Ds 


2,h.+ ee 
33”. ves 
1,4.+ i 
3,4 + ar 
2,n.+* ee 
I,n,+ Hest 
ayn, % toy 
2,n,t Weed 


Pe ee) © are 
35n. ; 


m. ) 
3,m, + j 
aaa ) 
2,m, ; 
cia. ) 
35m; ; 
Ine gigee Pane ) 
>) a j 


2,m,+ Way Zz) 
2,m,+ nM 2) 
1,m,+ ee 2) 
3,m;+ we ) 


yes ) 
2,m, é 


J 


ml (hat 1)y-A(L + L 
af 


+ 2 ) 
oon. + Oe 2a _ 
al 
J 


+ 


Vau(h + 


en 4 ) 
Poder @ ae n.,m,+ 3 
Ae alae 1) -A(L «+ ae yn, j 
J 


a(n s+ 1) =Z] 


: £ 
B(2 .+ 


Lh .+ 
A(ZL,+ : 


pate 


i) 


. L + 
B(4 + j 


2,n,+ nee 
2,n,+ a 
3,n5* ee 
1,n,+ eas 
Sy Use fey 


+ rt 
2,h, F 


Mo ee } 
2,m,+ j 


m. ) 
3,m,+ j 
Boise ts ) 
3 1 j 
om. ) 
3,m, : 
+m. ) 
I,m, + 





+ 


+ 


Uxe(n + 1)-z] 
Lah + 1)-z) 
ao 
a 
a. (Ue ny an Me 
Cho 1): 


len eee 
a : 


te bpua Die tae ae 
A(n.+ il 2,84 "5 ’ al 5 


Sete feasts iets 
"A(n,+ 4.4 3,h5+ 0, iv "5 


2 hee ea. 
A(n,+ oy 2,44 Ge 1 at 5 


| eee: a2 
A(n,+ se "5 dees j 


opi, ee ve 
A(n,+ hee ars a ; j 


1) 
ee he te ae ae 
A(n. + et 2,h.+ ns hs j 


‘B(n, + tt 
*B(n,+ a 
‘B(n, + a 
"B(n,+ oan 
*.A(L + n+ 

“A(L + ea 

"A(L «+ ig 


. ne 
A(ZL «+ j 


+ 
A(z 5+ n. 


ae. + 
A(h+ n, 


2 
+ 2,.m.+ m.+ 
A(L + ea oe he 5 


.>m.+ m.+ 2 
A(L.+ n+ 2,n.+ ae : j 


2,4 + tee 
3,4 .+ oe 
1 se 
3,4, +4 ses 
2,n,+ ei 


+ £ i+ 
1,n, j 


Jes ae 
3ini+ f 


+ £ .+ 
2,n. j 


sua! 1p eae 
3,ms+ m. 


i. + m. 


ae 
2,7. j 


) 


ak rice ba 


, ome 
2 oie 
L. 
ane 
} L.A 
ane J 


Lot 
2 as J 


Meets 
3,m, 

.+ ™ 
2,m, 

m.+ Mm 
3,m; 


1) Pia eo? 
i . 





and 
= Peet 2. fret + ooein, + m+ 2) 
a j a j a j 
+ A(Z.+ £.+ 2,n.+ n.+ 2,m.+ m.+ 2) 
1 J 1 J 1 J 
+ { Penman. © 26 ren 2 omer Ts) 
~ 1 J 1 J 1 J 


+ A(£ + a 2,n,+ = 2,m.+ ae 2) } (As) 


Note that if the expression ( 54) is used to evaluate the A's 
and B's that each of the terms in (A-4) and (A-5) above must be 


— Xv +U 
divided by a” * : 


2. <SI(H - E)'x? and <c| (H - E)IX? 

Each of these matrix elements is made up of 26 terms which in 
turn are made up of an integral of the form (53 a-d) where A's and 
B's are as defined in ( 44) to ( 47) and immediately following. 

Now letting X in (A-3) above be characterized by n,£,m, <S,C]|(H-E) lx > 
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where it 1S assumed Bens B. are chosen when the left hand side is 


<S| (H-E)|x> and A., B. when it is <C| (H-E)|y>. 


3. Determination of the Required Integrals 


Clearly not every term in the above matrix elements requires a 


different integral. 


Many of them are identical or have zero 


8O 





coefficient. For example, for the 34 element basis set (M = 5), 
while there are 595 independent elements each in H and S, repre- 
senting 52 terms and 4 terms each, respectively, or a total of 

more than 33,000 terms, only about 700, or slightly more than two 
percent, use different integrals and all 700 different integrals 
can be generated using the recursion relations proved in Appendix B 
by actually carrying out the evaluation of about 100 integrals for 
which f = 2. Similarly, in evaluating <S| (H-E)ly >and <C|H-E)|y>, 
where there are 34 elements in each expression, each element in turn 
made up of 26 terms, for a total of 1768 integrals required, only 
268, or slightly less than one sixth are distinct and again, using 
the recursion relations only 103 terms for wp = 2 need be evaluated 
explicitly and they can all be found from 36 different numerical 
integrals. Nevertheless it is the numerical evaluation of these 
integrals which overwhelmingly dominates the time required to solve 
the complete problem and hence the requirement only to evaluate 
those terms needed. 

TO insure that no terms were missed in computing the integral 
tables, preliminary computer programs were written to examine each 
of the matrix elements and to display in tabular form which integrals 
were needed to construct the matrix elements for a given basis set 
Size. Tabulated were all distinct integrals required for which the 
coefficient was non-zero in at least one case plus those required 
in order that the recursion relations could be used to generate 
the terms required for fp = 2. From this information the integral 
table programs were constructed so that only those terms required 


were actually calculated, and none that would be required was omitted. 
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APPENDIX B 
Integral Recursion Relations 


i. Recursion Relations fer uw > 2{ 51] 


Consider the integrals 


re ee Vie Daa? “ee 
A(V,A ,u)= Var ar jr¥ I, 12 F(r,x.) (B-1) 


and 


Ve2 A-2 U-2_ 
1 12 (B-2) 


B(V,A ,uU)= \ ar ar x 


Mo rewrite (B-1) in the form 


A(v,A ,u)= 


ae a 2 
] 2 ge 


2 = = 
+ - 2 r osQ ) d 
eee ey 5) 


using the law of cosines. This transforms immediately to 


POY AG) = A(Vt2,A ,U-2)+ A(V,At2,U-2)- 2B(vt1,A+1,uU-2). (B-3) 
Similarly B(v,A,u) becomes 


V-1 A-lu-4 at awe 
7 ry 12 yar ayaa 


2 
r,)cos8,, 19% 2 


te 


B(v,A ,u)= Bew+2 ,AWwe2)+ B(V,A+2 322 )- 2\x 1” 


Let the integral in the third term be designated I, then 


~1 A-1l_u-4 


= Vv 
ee NN tI 2 ) \ry r5 12 


2 Sas: 
29 : 
F(x,,r%,)sin 120% 95 


Now recal] that the angular part of the dr, dr., integration is 


simply sin®__d0 whence the angular part if the second term in I is 


y bevcrmes a [eee 
Ib 4 5S 
Sin 9_.,_d@O 


2 2 
_ = 2 
Ves I, 2h 100589 1.4) 1299) 





wee 
ae | eel ba Doe 2 “Z |sin’ 
a 9 
u-2 Flr, \ de, 5 nae ey Soa ae 
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which may be integrated once by parts to give 


2 4 U=-2 : 
- [ose 45 cos@, sin, ,do 


2 
ee . 


and therefore 


I = A(vtl,A+1,u-2)+ < B(v A, ) 


from which it follows immediately that 





w=-2 
BE); u)="— {B(v+2,d ,U-2)+B(v,A+2,W-2) -2A(v+1,A41,u-2)} ( B-4) 
2. Recursion Relations for u = 1. 


When uw = 1, (B-1) becomes 


V2 A-2 


A(V,A,1)= \ry a 


1 =n =s 
I) a hte 6 be 


F(x 
12 te 


if 


Now recalling that in the coordinate system in use ty Ee 3 SO 


ac 


1/r,, can be expanded in a series of Legendre Polynomials 





and the angular integration becomes 





£ 1 

= 6 T Ss I, 1 
_ wes \ P,(¢os@,,)sin9,,d9, , = = Gay P,(u)du. 

ry O ty aul 


But Pte) = 1 so the integral on the right can be rewritten 


1 
. Po (e)P, (4) du 


and by virture of the orthogonality of the Legendre polynomials, 


all terms in this expansion vanish except for £ = 0. Therefore 


Megan) = _ACMIO. 2) . (B-5) 
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In a similar expansion in the case of B(v,A,1) all terms but 


£ = 1 vanish, giving 
ace 
B(v,A,1) = 3 A(v-2,A+1,2). (B-6) 
3. Recursion Relations for ae a (VAY) and B. a (YAH) 
ce ee SE en hee 


Referring to (53a) 


Ag(VoA st) = AL (Vs HW) + AS (VAL) 
and applying (B-3) above 


A_(V,A,u)= A, (A+2,V,U=2)+ A, (A,v-2,4-2)- 2B, (A+1,v+1,u-2) 


+ Aj(vA+2,U-2)+ A, (V2 ,A 2) - 2B, (v+1 oe lem =2 ) 


A(V,A+2,u-2) + Ao (v+2,A ,lb-2) - 2B. (v+l A +1 ,U-2) 


and similarly with B_(v,A ,U). Thus A. and a satisfy the 


’ ’ 


normal recursion relationships for | 7 2. For »~ = 1 


BR (von .1) >= A, (A,v,1)+ A,(v,A 51) 


alt 


~ A, (A-1,V,2)+ A,(v-1,A .2) 
but 


A,(v-1,A a) = Be (vor x2) = A, (A 2) 


SO 


A, .(VAs1)= AL ((v-1,A,2)+ A 


. | 1 3A-lv.2)- A 


Th ee (B-7) 


where AJ 1S chosen when A. 28 tO bemevaluated and A. is chosen when 


A, is to be evaluated. For B. o(voA +1) the same procedure yields 
b) 


gd 
B o(voA51)= 3 {A 


oe (v-2,A41,2)+ Ay 3(A-2,v41,2) - A (A+1,v-2,2)} 


Ss aa @: reo 
(B-8) 


which completes the derivation of all the relations needed to 


evaluate the matrix elements listed in Appendix A. 
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APPENDIX C 
Rayleigh-Ritz Variational Method 


Let an arbitrary quadratically integrable wave function Y be 
expanded in terms of the eigenfunctions of a bound system 
YS) aca % 
i Xi (C-1) 
where 
(H-E)|X,> = 0 
then on the function Y the expectation value of H is 


<¥ | H|¥> 


X- 
> eee PN 
Eg May eliiine 


Zz 
Ela | oe (E=2) 


ig Eo is the ground state of the system, then E s EB. sO 
| 2 


Blige | es eel Ee 
O17 BN i a at 


whence 


assuming 


< ¥|y > 


It 
i 


But 


< ¥july > 


Ey 


thesexmectation of H on™’, so if ¥Y is found so as to minimize Ey 


then Ey 2 E, and Ey provides an upper bound on the value of Eo° 
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In general 


Seas ¥|uly >/<y|y > : (C-3) 


If now a trial wave function an is chosen and expanded on a suitable 





basis 7. 
* N 
ee i ee (a) 
1=] 
whence = | | - 
Debes H b.b 
: J is es 2 a 5°45 
E = anh, ee eae = 157) eee 
16 
ZS b.b.<n.|yn. > Y  b.b.S 
i,j 23 j i. es 
where Bee = <n lala? and S;. =<4,1n > Now choosing the b's 
to satisfy 
OF. : _— 
bbe = OF. Bae LS »N 
al 
gives 
OE, c » x 0) * 
BD. SR Stat E. dp7 % _bb.S,. = sh ; aE oelleeer 
; i,j 25 45 by eens ed 3 i,j 2 5 19 
from which arise the N linear equations 
7 ae eee ES oe - 
NE 5 Se, A ; 5 1 ,N (C=-5) 


which 1s just the ret part of the Harris expansion technique, 
wherein the smallest eigenvalue EY gives an upper bound on E 

The above method extends readily to the excited states [52]. 
This 1S achieved by choosing a trial wave function which is 


orthogonal to all the states lying below the desired one. Such a 


function must take the form 


N 
Yao -F X_<x_ lo > (C-6) 
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where ® 1S arbitrary and the ee are the lower lying eigenstates. 


Then expanding Y according to (C-1) it is Seen that the a. are 


zero for i less than N, the index of the state desired, from which 


mt tOlliows that E, = E 


Hylleraas 
automatically 
being greater 


oyocem. Thus 


ay 

and Undheim Cea have shown that this condition is 
tI : : 

met by (C-5), then : eigenvalue of that equation 


than or equal to the energy of the al level of the 


it 1S seen that an important bonus in the Harris 


expansion, method is a cataloging of the energy levels of the bound 


system formed by the scattering atom and the incident electron. 


Included in this cataloging are the autoionizing levels of the 


bound system which give rise to many of the observed scattering 


resonances C44] . 


87 








APPENDIX D 
Discussion of Numerical Calculations 


1. The Eigenvalue Problem 
em ee eee | 
From Section II C.1, a solution for eigenvectors and eigen- 


values 1S sought for the equation 


(H_ - AS)e = 0 (D-1) 


Since S is not the identity matrix (D-1) is not in the form usually 
associated with eigenvalue problems occurring in physics where the 
Space is orthonormal and S reduces to the identity matrix. This 
generalized eigenvalue problem has been extensively discussed in 
the mathematical literature (eee |. In the present case, however, 
11 1S possible to reduce the problem to standard form and treat 
it by the simpler methods available. 

ope ter first etae Matrix 5, Mts iclements are-the inner pro- 


ducts of the chosen basis vectors 


SSeS eee ; D-2 
$15 7 <X,IX, (De) 
Since the basis functions are real Sas = oe and therefore there 
exists a unitary transformation which diagonalizes S [53], Physi- 


cally this 1S equivalent to transforming to a new, orthogonal set 


~ 


of basis vectors. In this new set 


| = aut 1> 
oe Ge xix! 
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and Si 5° ORfor sal) wee But diagonalizing a matrix in this manner 
means that the new matrix has as its diagonal elements the eigen- 
mmrues Of S. Therefore it can be concluded that all the eigen- 
meas Of S are positive and hence S is positive definite. 

Now consider the positive definite quadratic form associated 
Meet S, X °S-x. There exists a real non-singular transformation 


Meewcea VOetor y Stich that x = U -y and such that (Ref. 53, p. 337) 


i -T -] _ —“ =9 
Ve SU Sy *y (D-3) 
calitea@ the "“cannonical form" of x +S *x, where gq” = @ >)". Hence 
2 =i) 
tee Solna ye (D-4) 
Zeeee 1 is the identity matrix, So 
See (D=5) 


In other words, a positive definite symmetric matrix can always be 
factored into the product of a real non-singular matrix and its 
Pecrspose, But because S 1S Symmetric, this factorization is not 


wieque, L£6r 


represents n? equations, only n(n+l)/2 of which are independent. 


Therefore one is free to choose n(n-1)/2 of the Uy Ambir ar ti yi 





5 


Unless one of the new basis vectors is null, which would mean 
the dimensionality of the space is reduced - a possibility which 
can be excluded from the present problem. 
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and it. 1S convenient to choose 
— ae 0) £>m 


for then U becomes an upper triangular matrix (one in which all 
elements below the principal diagonal are zero). 
With the above choice of factorization of S, (D-1) can be 


rewritten 


ue x0 (D-6 ) 
eee ae ae ae =r 
where x = U-c. Since H 1S symmetric it 1S trivial to show U ~-H:U 
1s also symmetric, for 
-T -1 -T -] -] -T 
= ; = 2% tl 
= so Dae ike ik a3 1k i,k im ki jk 
=p -1 -T -l 
= >: ; <— e . 
i,k jk 5 Fem oe ae 


and the problem has been reduced to the standard problem of a real 
Symmetric matrix in an orthonormal space. 

The advantage of the form of the factorization of S chosen is 
Eiatepecause U 1S triangular its inverse 1s particularly easy to 
obtain. The factorization 1S done by the square root method of 
Cholesky [54] and ut is then found by solving the algebraic 


equations [55] 


yy, = 
k ay = 1 ee 
>} UlU..= 0, i<k (D-7) 
j=l i) jk 
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Because the range of values of the elements of S and H was 
often several orders of magnitude (sometimes as many as fifteen 
or more) an additional check on the accuracy of ow as found by 
the above method was used. If the ch found above is taken as a 
first approximation to the exact inverse, then the following 


iteration scheme can be used to improve the value of the inverse [56]. 


mee ; : : = 
ext Sc; be the 1 0 approximation to the exact inverse, U - and 
@etine the error matrix 
Re f= G.-ue (D-8) 
= = a Fo 
Then 
; — i 2 . ¥ = Gs = : ° ’ 
and 
; = = at a Cs 
pe = Eee (l= coauiec.) *u 
2 eZ p=9 
= (I - C.:U) =R. (D=9) 
= ke 2 
x eet 
= (R°) 
where 
ee ee Oe 0 


and om 1s taken to be the value of ee found by the Cholesky 

wethod above. Clearly, if the norm of ies is less than one the 

method converges quite rapidly, and since the value of Se is 

already quite accurate, WR II<< ies In practice no more than one 

iteration of Cc. waS ever required to make the elements of R< 10" oes 
The eigenvalue problem was solved uSing the Jacobi variable 


threshold method [41,55]. This method proved to give the most 


accurate results although it is known to be considerably slower 


pl 








than other methods tried. The accuracy of the eigenvalues was 


verified by the following method. For a given eigenvalue, A the 


(ie 
determinant |H-\ S| was formed and compared with the determinant 
when aT was replaced by aate® where 6 was typically 1077 ( compared 


to dey which was on the order of one). Then Srp was considered cor- 


rect to six places if the respective determinants were related by 


Ju-(A,, + 6)si<|H - A si<|H - A, + 5)sI (D-10) 


and the first and last terms above were of opposite sign. Accuracy 
of the elgenvectors was not so clearly established, though a 
Similar procedure to the above was followed. The product (H-AS)c 
was formed and it was noted that as A was perturbed as above, all 
elements in the resultant vector changed sign. While this indicated 
that the eigenvectors probably are reasonably accurate no quanti- 
tative bound could be piaced on their’ error. Further indication of 
the accuracy of the eigenvectors can be inferred from the close 
agreement of the phase shift results for hydrogen with the accepted 
values found by Schwartz G12 
2. Numerical Intearation 

From the beginning it was felt that the simplest and most 
reliably accurate means of evaluating the necessary integrals 
numerically would be some form of Gaussian quadrature C57]. Any 
numerical integration scheme attempts to replace the integration 


by a finite sum so chosen as to make the error tolerably small 


b nN 
( i x) cr > w.2(x;) (Daia5) 
ah 


a2 





where the w. are weights to be assigned to the function evaluated 
at the sampling points Xo in other words a weighted average of the 
function at the various sample points. If the Sy SUES fixed, equally 
spaced points then the n+ 1 Ww, can be chosen to define a polynomial 
Sreor@er n wisth which to approximate f(x), and if f(x) is itself a 
polynomial of order n or less, the integration will be exact. If 
however, one is free to vary the Xs also, there will now be 2n + 2 
parameters available and it will be possible to define a polynomial 
Srecearee 2n + 1 wath wvhich to approximate f(x). Thus, given the 
values of x. and Ww. in each case, the same amount of labor results 
in a more accurate integration in the second case. 

Adopting this approach, one can approximate f(x) with a 


Lagrangian interpolating polynomial of degree (n + 1) of the form 


nN 
‘ae ca : : D-12 
£(x) = o(x)F(x)= E_o Ly (x)a(x)F(x,)* Ry (x) ae 
where 
Nn 
ee) See (D-13) 
. j=O (x.-x.) 
jAi 
and Rf) is the remainder term, given by 
nN 
R(x) = TU (xx, )P 12) (2) g(x) /(nt1)!, a <B> bd (D-14) 
i=0 


The reason for factoring f(x) into g(x)F(x) will become apparent 


e a E + 
below. Now 1f F(x) 1S a polynomial] of degree 2n + 1 then pi? > 


~~, 
uy 
ad 


must be a polynomial of degree n. Letting 


(n+l) 
ees ) = gq (x) 


(yeeie)e: n 





where q_ (x) is a polynomial of degree n, (D-14) becomes 
n 
RCs) = a. (3) TT (x-x,) 9(x) (Dt) 
1=0 


Integrating (D-12) gives 


b n Ge b - 
J 9) F(x) 0% peg ete, cx + \ 9 x) q(x) IT oh (D-16) 


Now (D-16) is of the form (D-1l1) where 





b (Sic 5) 
ee Vo OTT 5 dx (D-17) 
a i 
#8 
except that the Xs have not been chosen. The object is to so 


choose the x, that the second term of (D-16) vanishes. This is 
easily done if q, 6%) and IT (x - x5) are «panded in a series of 
polynomials orthogonal on [a,b] with respect to the weight 
function g(x). A weight function appropriate to a particular 
set of orthogonal polynomials is often a natural factor in the 
function to be integrated and therefore the particular expansion 
to be chosen may be dictated by the problem, hence the utility in 
carrying the factor g(x) through the derivation. 

Let the polynomials chosen for the expansion be designated 


P(X) whereupon 


Qe 
~18 
TE (x25) = Edgy (2) ae 
i=0 
and 
n 
Bie) > ze Sie) (D-19) 
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Combining these two expansions in (D-15) and applying the ortho- 
gonality condition 


b 
\ o(ry GaP, (s)dx = 0, £ 4m (D-20) 
a 


reduces the error term to 


‘ : : " otxp CP at)? 
590909 TT boxy )ax = Z_ bse, 9%) P(x) ax ~ 


Now if the x, are chosen to be the zeroes of ee 16%) > then the 


+ 


expansion (D-18) reduces to 


n 
ROSS) se eae 
1=0 
Te .:, De =O OSjSn, whence (D-21) vanishes, and the integration 


1s exact. Thus evaluating the function at n+l carefully chosen 
points suffices to integrate it exactly if it can be expressed 
in terms of a polynomial of degree 2n + 1 or less, a considerable 
improvement over the equally spaced interval methods, as was 
promised above. The one remaining problem is, then, to evaluate 
tive Wis and find the x, "Ss. Fortunately, this has been done for a 
large class of orthogonal polynomials by Stroud and Secrest [58]. 
Although the integration limits and the form of the integrands 
(49), (51), (56) and (57) suggest the use of the Laguerre poly- 
nomials for the numerical integrations, the almost periodic nature 
of the functions and the ever increasing amplitude in the absence 
of the e factor dictated the use of a finite interval formula 
where several intervals could be integrated separately and then 
added together to give the final result. Since the factor oa 


forces the integral to zero at large x, examination of the integrands 
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indicated that an upper limit on the integration of x = 100 would 
introduce a negligible error from the neglected region in the worst 
possible case. Because the variable change to conform to the con- 
vergence limits of the polynomials used in the expansions (D-18,19) 
a Simple and because the weight function is unity the decision 

was made to use Legendre polynomials in the integration. 

Because of the almost periodic nature of the integrands, when 
the instantaneous frequency of the sinusoidal variation is large the 
areas under the integrand above and below the axis are nearly equal 
and substantial loss of significance occurs when these two values 
are subtracted. Often aS many as eight to twelve significant 
figures can be lost due to subtraction, hence great care must be 
taken to avoid the problem. The approach adopted was to seek a 
function whose exact integral 1s known and which approximates the 
desired integral closely enough so that the difference between them 
is small, then evaluating this small difference numerically and 
adding the result to the known value of the approximating function. 
This method can be compared to the measurement of the difference 
between two large quantities, say frequencies. Often the difference 
can be determined quite accurately even though the absolute magni- 
tude 1S known only approximately. 

In executing this approach to the integration of the functions 
at hand the first step 1s to find all the zeroes of the integrand 
between arbitrarily chosen minimum and maximum values of x. These 
values are then used as the end points of the subintervals in an 


integration of the form (D=22) 


cl A ae Zs 4 100 
d=) £(x) dx+ iad, LCs) ng) Jase) ns Cas), f(x) dx 


n+] 
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where the z, are (ie@n + lezeroesso1 1(x) found above and the h. (x) 
are functions whose exact integral is known and which approximate 
f(x) over the interval indicated. The h, (x) are chosen as indicated 
in Figure 23. Note that h. (x) 1s amplitude modulated from segment 
to segment so that f(x) - h. (x) 1s always mostly positive and 
smaller than f(x) by at least an order of magnitude regardless of 
eimeesign Of f(x). It is this difference function which is inte- 
grated numerically, and since h. (x) isa Gitierent tunetaonm 1 or 
each half cycle of f(x) the numerical integration is done in 
segments between zeroes of the function uSing a 32-point Gauss- 
Legendre quadrature on each segment. 


The form chosen for h. (x) was 


(D-23) 


where p is +l if f(x) 1S negative and -l if f(x) 1S positive, 
emit. Or 1. is to be found and §5= O72(7210) if I. or 1, is 
Ss e 7 3 


to be found, and 











TT 1-P 1-P | ) 
= = Paden + ———)z. D=2 
A 5 seemed al (1+ 5 2544 2 5 )z ( 4) 
pH Sle 
so chosen that sin(Ax + B) vanishes at the end points of each 
segment and nowhere in between. 
When chosen in this manner the integral of h, (x) alone 1S 
exact and has a rather simple form. If (1 - ye is expanded 
fomr in@egrals of the form 
Pt l -a,X 
Sy =| sin(A,+ Bye x’ dx, £ = 1,2,3,4) (D-25) 
cae 
1 


OF 








Figure 23, 


Example of Method of Choosing h. (x) 


mesuity where a, is 1, 1 + B, 1 + 28 and 1 + 38 in turn. These 


L 


integrate immediately to 


yr “A,Z. : =a P 
Sp, = Yd oo e ‘ as aimless ' a 
Jal 


yf D-26 
(-0) J(y+1-3)! sei : — 
where 9 and 9 satisfy the relationship 
My 
~a + (-1)°A = p(cos® + (-1) “sin@) 
and finally 
7+) 
\ h.(x)dx = P(1+=-)(S -38.+ 35.- § ) (D-27) 
; i LO eed See: 
i 


98 





Expressions (D-26,27) are valid for y20. During the course of 
evaluating the phase shifts, however,it is necessary to evaluate 
the various integrals for Y = -l. For this case a slightly dif- 
ferent scheme must be used to evaluate h, (x). If one replaces 


fel = B75 in (D-23) by the integral identity 


_ PX tone 
Ue ae = B ( e oe (D-28) 
Be 
O 
then 
z : ; -x =Bso 2 Bxy (D-29) 
h.(x) = B(1 + zz) \ sin(axtB)e “(1-e '*)“e Fay 
1 10 
O 
Carrying out the x integration first and making the appropriate 
variable changes yields 
Sry) i Mh Ba ill 
S, = SS eee ae (D-30) 
£ 2 Z 
BL+1 w +A 
and finally 
all p 
ies alse = =: Te (lier oe) eal ie 2 Ga S 
\ heals ( i0? ( O 1 2) (D-31) 


a 

Unfortunately (D-30) must be integrated numerically also, however 
1t 1S a relatively slowly varying function of wand is always 
positive so integration by a 32-point Gauss-Legendre quadrature 
gives more than adequate accuracy. 

To complete the integration of (D-22) the first and last terms 
were evaluated using Single 512-point Gauss-Legendre quadratures. 

Sumpeesthe- accuracy of the numerical integration is crucial to 
the accuracy of the phase shifts calculated for helium ion scat- 
tering, a great deal of effort was expended to insure that the 


methods were both correct and internally accurate. 
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The first and most obvious check is to increase the number of 
sample points in the numerical integrations. This was done by 
increasing the quadrature method on each interval to 64 points. 
Changes to the integrals on doing this were typically in the ia 
to a digit, indicating that the differences were principally 
due to round-off error since the maximum number of digits available 
was 16, 

The integration scheme described above grew out of a similar 
method which was used extensively prior to the present one for 
the same calculations. The present method was adopted in an effort 
to extend the useable range of the parameters involved in the 
integrands. The problems with the former method arose principally 
because it was not sufficiently accurate when the oscillation fre- 
quency of the sinusoidal function was high. However in the region 
where both schemes were adequate (which proved to be almost every- 
where) the integrals were consistently reproduceable to 12 to 14 
places. The choice of integration intervals and of the h(x) me) 
the earlier method was sufficently different that 1t constitutes 
an almost independent check on the internal accuracy of the 
integrals. 

Three different methods were used to make the integrals exactly 
integrable so that the numerical results could be compared directly 
with known answers and to provide a check on the correctness of the 
Cypemessivons. The first was to let k ~ © while holding the ratio 
Woe ceovistamt. This forces the logarithmic terms to zero and re- 


duces the integral to a standard form which can be integrated 


directly. The direct evaluation was done on a Hewlett-Packard 
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Model 9810A computer which gave results to ten significant figures. 
These results were then compared with the numerical results and 
except in the cases where k/@ was at its extreme limits, agreement 
was obtained to eight to ten places. This method was tried only 
with the older integration scheme and so some decay at the limits 
of k/@ was expected. In addition, the transition from the form 
of the integral when k was finite to its form as k ~ © appeared 
smooth and with no obvious discontinuities or other strange behavior 
in the transition region. 

A second check took advantage of the fact the computer program 
was written to allow for variation of the nuclear charge, Z. 
Although the intent was to limit the value of Z to integers, 
nothing in the program or the mathematics forbade varying the 
charge continuously from one to two. When this was done, it 


was again seen that the result for Z = 2 transitioned smoothly 


I 


maco that for Z 1 with no discontinuities as the charge was 


varied. 


I 
—_— 


Since for Z the integrals all become exactly integrable, 
the computer program did not do the calculations required for the 
numerical integration in the case of Z = 1, but skipped to a sub- 
routine which evaluated directly the exact integrals. The final 
comparison method was to circumvent the exact integration routine 
for Z = 1 and do these integrals numerically, using the method 
described at length above. Again the agreement between the two 


calculations was excellent, the integrals and the phase shifts 


typically being identical in the first twelve to fifteen places. 
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MWe final check was to test the sensitivity of the entire 
Saveculatyvon to perturbations in k. It was noted that even at the 
extreme limits of k/Q@ where accuracy was expected to be poorest, 
only some of the integrals showed evidence of loss of significance 
from the subtraction problem mentioned above. For these integrals, 
however, a small perturbation in the value of k may have a rela- 
tively laroe effect on the integration and if the inaccurate 
integrals dominated the computation then they could significantly 
distort the results. As before when this perturbation was arti- 
ficially fed into the problem its effect was well within the limits 
considered acceptable (results unchanged to six places or more for 
a perturbation of k in the seventh place) except at the extremes. 

Since none of these tests indicate any serious problems with 
accuracy it is felt that the errors inherent in numerical inte- 
gration have been kept well within tolerable bounds in this problem 
and the results can be relied upon with considerable confidence. 

It also proved to be true that those regions where accuracy was 
lost were easily avoidable. 


pre evaluation of IT (1 = ia) 


From the definition of the gamma function 


r(P) = \ ie ea (D-32 ) 
O 
one finds 
— . 
Y(l-ia) = \ ex Fax (Dae) 
O 
but es x, so (D-33) can be rewritten 
"(l-ia) = \ e *cos(a&mx)dx-i \ e “sin(alnx)dx (D-34) 
O O 
sags te eg Sale: 
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whence 


argue 1 -iayie= ects 7 (D-35) 


and the evaluation of the complex gamma function is reduced to 
evaluating two real integrals. The integrals have to be evaluated 
numerically, however, and are not suitable for this in the form 
(D-34) because the rapid oscillations of the integrand near the 
Origin could contribute to a substantial error in the numerical 
integration. To overcome this problem note that an integration 
by parts of each of the integrals (D-34) introduces a factor x 
into the integrand and this has the desirable effect of reducing 
the contribution of the integrand near the origin. Each re- 
petition of the partial integration introduces an additional 
factor of x into the integrand. Carrying out this process three 


times transforms the integrals (D-34), after considerable algebra, 








into 

Be Sth ena ee ae | x30 {8 (2-1) sin(atx) -(a7-11) cos (ax) fax 
(1t+a™) (4t+a™)(9+a™)"O 

i, (D- 30) 

n 

A = —~—— AG 36 aoe {(a?-11) sin(ainx) + °(a7-1) cos (a’-x) fax 

(l+a 2) (ya?) ) (9+a~)MO 
(D-37} 


which forms are now quite suitable for numerical integration. 

The integrations of these functions were carried out uSing a 
256-point Gauss-Legendre quadrature over 15 arbitrarily chosen 
intervals. To check accuracy the results were compared with 
another program which calculated the complex gamma function using 


a Pade-power approximation method and claimed 9-place accuracy. 


103 





The internal accuracy of the integration scheme was checked by 
comparing the 256-point integrations to 512-point integrations 

on the same intervals. This comparison indicated accuracy to 
about 13 places and these results were consistent with those found 
by the comparison program. In addition the numerical integrations 
seemed to take somewhat less time than was required by the power 
series method and was not subject to slow convergence as the ima- 


GQinary part of the argument became large. 
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COMPUTER PROGRAM 


PROGRAM PHASE MK IV SEP URAC Ta 
MABE ICIT REAL*®8 (A=Hy0-Z) 
PRAEGER RO 


EXTERNAL CGARS,FACT,DCOS,DSIN,FON,FCN 

Seep oomMP LE POINTS AND aE IGHTS FOR 256 AND 512 POINT 
LEGERDRE INTEGRATIONS. 

COMmmON7 tL 256/X256(128) ,A256(128) 

GO, 691274512 (256),A512(256} 


PreetORATION INTERVALS FOR GAMMA FUNCTION INTEGRATIONS 
CUPmUN7CO/DPC{( 15) ,0NC{(15) 


MPISCELLAMEQUS CONSTANTS UScD IW VARIOUS SUBROUTINES 
Cmermen7, CUONST/SPRIVP AZ. FRLsFRZ,ASCL,P1,CA,sCT AA 


a ALPHA AND AND ETA TO THE INTEGRATION AND MATRIX ROU- 


N 


COMMONS TR4ASAL-A,ETA 


Peeeoes K TD THE WRTSGRATI GIP ROUTINE 
CUMMON/ WAV PZWAVN 


Pes seS WUCLEAR CHARGE 
COMMON ATNO/Z 


Saad VARIAEGES. “T4t41S COMMON USED IN INTEGRATIONS FOR 
RHG = - l 


COMMON/MCNE/D1 02,03 


USED Uehiro (He LIMITS OF THE RANGE OF THE INTEGRAND 
@VER Whee [0 FIND THE ZEROES 
COMMON/INT/XS“M4FF, BU sXLM 


ea pee eS THIS COMMON USED IN INTEGRATIONS FOR 
> - 


COMMON GEZRQ/94,05,06,110,120 
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merCES FOR MATRIX ELEMENTS 
COMHON/ INDEX/NINJsLIstJsMI,>MJ,18 


meLUE OF EXPONENT OF X IN BOUND-FREE INTEGRALS 
COMMON /RHO/RO 


pero neGe FOR BOUND-FREE INTEGRALS 
DIMENSION AS(1452),885(1452),AC(1452),BC0(1452) 


pgp e FUR THE V4LJES JF THE EXPONENTS OF Rl; R2 AND R12 IN 
Mde HY LEERAAS FUNCTIONS. LOCATION OF THE REGICNS OF ALPHA 
Mem oe son PPPED IN THE CALCULATICN 

Oem NS TONeN( 95) +LI95)+M(957,ISTP(10) 


Pea stiveo'’s FOR EXACT ENERGY CALCULATIONS 
DRRENSTONTX(5) .Y(5) 
c DMMeENSTON STATEMENTS FOR PHASE SINGLET, N=7. 


STORAGE FOR BOUND-BOUND INTEGRALS 
Sareno LOW £15508) 5 615508) 


See ee ce FOR MATRIX ELEMENTS AND EIGENVECTOR ELEMENTS 
DiiiegS! ON H070.70) .-St7JIs,72),C070,70) 


STORAGE FOR EIGENVALUES 
DIMENSION WAV( 70),cNERG( 70) 


Piece FOR UPPER TRIANGULAR ELEMENTS OF H AND S$ AND FCR THe 
Pemeeeoe OF THE TRIANGULAR FACTOR OF S$ 
Cewe NSION TH(2485)51S(24685) 
E@UTVALBMCETHTI), ASTI) ), (H(2451),BS(1)) 
EWOUVeWenCat ol yee wires (51 2451),BC(1)) 
e 
Pewee FOR HE AYLSR&AS FUNCTIONS AND THE NUMERICAL IMTEG— 
RP PON eres STOMED IN A ScCOUENTIAL DATA Set 


READ(S8,1lOOOTN,L.4 

Rew Bl O02). DPC,OMC 

RESO (3,1 OM@d (X51 2011), A512(11)5 11-1, 2556) 
READ(8,1O001)(X256(KK),A256( KK) ,KK=1,128) 


Rievmeweocre OF DATA FUR PRESENT CALCULATION IS READ IN FROM 
PIMCHED CARDS SUBMITTED WITH THE HAIN PROGRAM. THIS DATA 
Realms) ae? FOR [4E eNTIRE PROGRAH: INDeX FOR NUMBER OF 
Viteeopat loans, TOTAL SPIN; BASIS SET SIZE DETERMINED BY THis 








Meee Xs NUMBER OF ENERGIES TO BE INVESTIGATED (SET THIS EQUAL 
Memeo PP lic ENTIRE ELASTIC RANGE IS TO BE COVERED}, 
meCLEAR CHARGE 

READ(5,1021) NSceTs,IS+sMATRIX,INO,Z 

Petes. 6On0) Gi» TO 88 
Meme seL SPIN IS GQNE (TRIPLET) SCREEN OUT THOSE FUNCTIONS 
Bere wEEDED FOR THE CVALCULATION 

K = O 


DO 109 KK=1,95 
IF(NC(KK) .EQ.L(KK)) GO TO 100 
K = K + 1 
NCK) = N(KK) 
L(K) = L(KK) 
M(K) = M(KK) 
Metim@cO.tATRIX) KK = 95 
1J0 CONTINUE 


See eeOUS COMSBTANTS WHICH DETERMINE NEEDED INTEGRALS 
So wmex = AeNSET + 4 
WAX = NSET + 5 
MAXI = NSET + 3 
Pier = MATRIX=PATRI X 
fl ZE =e xetATRAX + ae) 72 


mmo THE REQUIRED BOUND-BOUND INTEGRALS 
CALL TBLI(MAX+7A_B) 


SET CGMSTANTS N&€GDED FOR THE NUMERICAL INTEGRATIONS 

tale) woe 2 OOD Oo to) 

ASM = 0.50-2 

Fre U0.2)D-] 

BO = 0.1002 

XLM = 0.45502 
Peewee ON TA FOUR THe RANGE VGE "“AEPHALAND ENERGY 912 
TIGATeED. 

DO 139 IW=1.,INO 

READ(5,;1008) IALO, IAHI,1AINT »NUM,FCT,PLO,PHI 
VERSUESTOEV APPA TCO BE SKIPPED 

pe Oto y 1005) -1STP 
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Peeeoetomwetil | SeeRCt FOR CLOSEST EIGENVALUES IN EXACT ENERGY 
Memeo ies. IF ND EIGENVALUE IS FOUND WITHIN THIS EX- 
evel WANGE THE CALCULATION IS STOPPED FOR THAT ENERGY 

Pec Aw = MHI + 0O.1D-1 

Mig@ee.f = PLO - 0.1D-1 


ene Petre OoveRS DESTRED RANGE GF ALPHA FCR DESIGNATED 


1@2 OO 10) LA=1 ALO, TAHI,TAINT 


me tie Om VALUE OF DESIRE Oi S 4820, aide ENTER CLASTIC 
Meer ter ING KeGlON WILL Be INVESTIGATED. SKIP THE EXACT 
ENERGY CALCULATION 

fer Loe cee GOI) GO TQ 50 
Sem lide evel UES CF ALPHA AND K FOR THE FIVE POINTS GN THE 
eee rery Cleesocst TS THE DESIRED VALUE OF K 

Reape lUM6) LY (ide Xt1) + 1=1 959 

faa = PLO + U0 DS 6 
Use A CUBIC IRTERPOGLATING POLYNOMIAL, FIND AN ESTIMATE OF 
ferme en We OF PALPRA REQUIRED TC GIVE THE DESIRED VALUE OF K 

70 C@EL SP ths Y,5, XoewT ae YINT) 

ALFA = YINT 

GO TO 5l 
eee fo GEE. IF A SSG4ENT OF THE ALPHA RANGE IS TO BE SKIP-— 
eee em USED Where THE EXACT ENERGY CALCULATION IS USED 


50 CO 103 I=1,NUY%.2 
Viet itie eer e C1): TAS ISTP Cle 
103 COMA IMUE 
hea = FORe eee OAT 1 A) 


Eeeweeetc THE H AND S MATRIX ELEMENTS 
51 DO 16 J=H1-.fiATRIX 
ee ee / 2 
DO 10 {=l-sJ 
NI = NCI) 


NJ = NJ) 
ee = Lott) 
LJ = LJ) 
Al = #1 ) 
HJ = M(J) 
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Cree ELEM(HH, SS:A+BbsMAX) 
Pest + J 


H(Ie«J) = HH 
Mitts) = HH 
S@F,J) = SS 


het oti) = SS 
ie “CONT TINUE 


Mmemere sr ACTORING S, SCALE IT SO THE LARGEST ELEMENT IS < 10. 
Mm l1> tHe SCALE FACTOR 

Geet SCE(7TS. 1S 12ZE,F) 

WRITE(6,1006) F 

Cais SUBY(1S,Fs TSIZE) 
moon Ss 1NTO ITS TFRIANGULAR ELEMENTS AND INVERT THE UPPER 
(WremesGULAR FACTOR. RETURN IT IN TS 

Gee, INVERT(TS+:HsS+CesMATRIAs ISIZE) 


pee Ure (—-1 )SHEUS=(-1). RETURN IT IN TH 
Comer MmeRD (THs TSsC sMATRIXGF) 


CED SUC e Vieticom rules LN UPPER 
Veen Puts e ASc 
I 
if 


oe J 
DO 20 I=1,J 
[ieliesal, wet) o 
HT = H(1I.J) 
H(1I3J) Taki] 5 ) 
Hiss) =) Aer) 
TH( TJ) = HT 
Say? = 0.009 
20 Side) = TStId) 
VRACE = ODODO 
DO 11 T=leMATRIX 
MAG Re TRACE ttt J 5 T ) 
MRICS CS. 2013) TRACE 


FimiD EIGENVALUES AND EIGENVECTORS OF TRANSFORMED H MATRIX 
CALL OJACVT(H»sMATRIXs1 sENERG,C,MATRIX) 
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TRACE = 0.050 


FIND 
IOUS 


A Ui i ee E> AND COMPARE WITH TRACE OF H PREV- 


AY FOUN 

DO 17 T=1,MATRIX 
Piette CE > TRACE + ENERG(]) 
WRETE(Gs 1015) TRACE 


Peete weewoekvECTORS $O LARGEST ELEMENT 1S < 10 
She UYSCUCs LIME) 
weet oMry(C,yF, LI) 


mOLTAPLY THE ORTHOGINAL E 
Mie EYGENVECTORS OF THE O 
Geel -M@PROCS yy Hs MAT 
MRI EL64 3002) ALFA 

WRITE(643004) ENERG 


IK = O 


I BY = 1) eo: ND 
E L 


7 Oe 


Sees CAMs The ELGENVALUES TO FIND THOSE IN THE 
ED EWERGY RANG= 


DO 15 TE=1,.MATRIX 
WAV(IE) = 0.0900 
ct = ENERGCTE) 
BOOT =OeZDieeE. + 257 
NetR@@T .LE.0.JD0) ROGT = 0. 0D0 
PINT = DSQRT( ROOT) 
EROQPINT.SEE.PLO)-OR.(P UN WaT .Pai co 10 2] 


‘ GIVIWG NUMBER OF EIGENVECTORS IN THE DESIRED ENERGY 


IK = Ik + 1 

WAVCIK) = PINT 

mlees iit DESTRED RANGE RELOCATED "IN FIRST ADDRESSES OF 
MALGE VECTOR 

PMERGCIK) = E 


= 
a 


E 
EN 


Cites CORRESPIND I 
WATE COGMMN- OF E 


DO 12 K=1sMATRIX 
mee tiicpikelee= C(K, IE) 


N 
| 


QQ 


eas RELOCATED Si 


PO 





ION 
NE 
EF = 


—™ C4 


ASSUME 
COG NED- “Ff 


MATRIX 


ee VALUE WITHIN 


Memexeet SNERGY FEATURE IS NOT BEING USED SKIP THIS SECTION 


Crete ro. EQ. Je0D0) 


6) 


G2 


61 


63 
65 


GO TO 15 


EemeLUE BEING EXAMINED IS OUTSIDE THE EXTENDED RANGE 
Pomme NEXT ONE 

amet NT ole PLOEXT IOOR- (PINT .GIePHIEXT))GO TO 15 
Pe eeUe TSseCUrSIDE, THE DESIRED RANGE SUT 

Be JNE>, PRINT THE VALUE ANO REPLACE THE 

ee Mele SEGMENT CF THE TRAJECTORY UNDER E 

Peeve. THEN.Ge BACK AND RECALCULATE AN ] 

meeor THE ReE@e@IRED ALPHA 

WRITE(6+1014) IE,PINT 

KNUM = O 

oie INT .LTsXt1)) KNUM = 1 

IFCPINT.~GT-X(5)) KNUM = 5 

PeWeR NWO. 6Os1).0R. (KNUMLEQL5S)) GO TO 65 


Bis 60 T=2,5 
Meee tel t.Xl J.) ). ANOS. (PING. GI oxti-1))} 


CONTINUE 


Pet (KOM W.LT.1) -ORs(KNUM.GT.5)) 
DC 62 J=1,.4 


Lieo= 5 = 


I 


het Il. Geary) 


CONTINUE 
X(KNUM) = PINT 
Y(KNUM) = ALFA 
GO TO 79 

KNUM = KNUM — 
DUT 6S IT=zy5 
Pere ee. nN) 
PFC i .LE.KNUM) 
CONTINUE 
X(KNUM) = PINT 


pal Gl Uae ca | 
Y(IL+1 


1 


vt —)) 
YUl= 


) 
) 


X(IL) 
VOTE) 


X(1) 
alee 


CGOe TG 15 
COmTarol 


KNUM 





Y(KNUM) = ALFA’ 
GO TO 70 
15 CONTINUE 


ee ELGENVALUES WERE FOUND IN THE DESIRED RANGE, 


NEXT VALUE OF ALPHA 
Peis LOO) GO 10 101 
Were 1 EG, 3005) ALFA 
RITE(6,3004) (WAV( LP) »LP=1,I1K) 


FOURTH LOOP. bBeolN THe PHASE SHIFT CALCULATION 
DO 113 NN=1,I1K 
E = ENERG(NN) 
WAVN = WAV(NN) 


mest fF lieO THE COULOMB PHASE SHIFT. ETA 
AA = (Z - 0.191) /WAVN 
ee 4 
RATICNS 
IF fAA.EO.0.0D9) GO TO 301 


ASQ = AA#AA 


CA =ASQ =e TED Z 

CT = 0.6D1*(ASQ =O. LOL AA 
Gs INTZ55(C5ARG,GR) 

Cee = CT 


Cl] 2051102 —A50 
Cael -bieh2>6( Co ARGWwG61 ) 


6% = —GR 
ETA = DATAN2(GI1-GR) 
Ga5 10302 


301 ETA = 0.0D0 
WRITE(641320) MATRIX«MATRIX,ALFA,EsWAVN 
Ga. 10:-203 


392 WRITE(6+1004) MATRIXsMATRIXe¢ALFA, E+ WAVN,ETASZ 


FIND THE REQUIRED BOUND-FREE INTEGRALS 
Bw -Crmet TBL2{ IMAXaIWAX], AS, AC.BS,BC) 
STERM = 0.0D0 
CTERM = Q.O0DO 


a2 


CO LOsthe 


: Peo P Steers) AND GO ON TO ihe BGUND-FREE INTEGS 








FORM THE 
DO 
ve 
NJ 
LJ 
MJ 
Gar t 
HSC 
Sy bere a 
nec 
CekM 

300 CORT INUE 


300 JM=1,MATRIX 
C(JMNeNN) 

NC JM) 

L( UM) 

Pe ry) 


HS*CC 
STERM + HSC 
HC «CC 


Chekiesd Fee 


— 
-_— 


pvereentE™RHE PHASE SHIFT 


TOL TA = -STeERM/CTERM 
TOLTAK = TDLTA/WAVN 
DLTA = DATAN(TDLTA) 


Be (lS See 0) eo (0 91 
perre to. 101) 
Ce iO 92 
Ol WRITE(G.1012) 
SZ MR ITEISs 1007) 
WRITE€ 7511009) 
Pet Oet a, Gl«0.000) 


COT) ite 


Neeeide sean St SHIFT |S NEGA 
PRINT AN ADDITIONAL DATA C 
DLTAP DEA, Pl 
WRITE(7s,1101L) ALFA,WAVN+-DLTAP 
CONT INUE 

CONTINUE 

CONTINUE 

See 

FORMAT (4811.+/+4711) 

FORMAT (2030.16) 
PURrman (809. 55/4709. 3) 


TIVE, ARB 
ARD Wel ta 


2 
eel 
PeO 


1000 
LOO1 
1032 
1003 


Is 


BOUND-FREE MATRIX ELEMENTS 


PHASE(AS: AC, BS ,BC +MAX],HS+HC, E) 


ALFA, WEVN;ITIDtLTA, TOLTAK;DLITA 


ADD PI TC IT AND 
Beiierse (Silt | 


EGmmee ('—"',5xX, "SOLUTIONS FOR ALPHA =*,FIL1.6,/;'=-*) 








1004 FORMAT('1',/,'-',10X.' ELECTRON — HELIUM + SCATTERING 
1PHASE SHIFT FOR',I3," X ',12;" HAMILTONIAN MATRIXe 'q/¢ 
eae ALPHA =',F10.71/.'0'+24X%,"TOTAL ENERGYs E =", 
eee / et 0'.25X, "WAVE NUMBER, K ="'.D24.164/¢!O%15X,_' 
econ PHASE SHIFT, ETA =!,024.16+/e83'+20%> "NUCLEAR 
Sormince’, Z =", F9.5) 

1035 FORMAT(1018) 

iglo6 FERMATI™I,1OXe'F = ',D10.3,/,t=") 

ewe G@RMAT('O'.5X.* ALPHA =',F9.604X%5'K =' 5 D2401655X_! 
1TANGENT (PHASE SHIFT) =',D24.16;/;'0',65X, * TAN(DELTA) /K 
eet G7 0". 66X%. "PHASE SHIFT =". 024.16) 


1008 FORMAT (80m 218, 213,F10.355X,2F 15.9) 

er ear Te «fe? —" LOX," TRIPLET PHASE SHIFTS.* } 

eee er eeAT(O'«/+'-',.lOXs'SINGLET PHASE SHIFTS.'*) 

ero Urea tT( LOA, "TRACE OF (Us*(-7T))*He (USF (-1)) =*, D24e16) 
Meera TOMMAT( '—" ,LOX.'*KEC' 12,5") =',024.167/,'-') 

mer> FORMAT(22X,'SUM OF EIGENVALUES =',D024.16;/) 

1016 FORMAT(F11.7,F 12.7) 

moi eer (*—".10X,15,.' CARDS PUNCHED FOR THIS JOB.*) 


moro eewwe: | Pe / 5 '=',10X,' ELECTRON — HYDROGEN SCATTERING 
Memeot SHIFT FOR',13-' X 's]2+* HAMILTONIAN MATRIX. aes 
Cee ve Oe PHA. = hE 965/53 ' 0's 24Xy' TOTAL ENERGY.) Eo = 7% 
SU25516+/,.°'O"',25Xe "WAVE NUMBER 9 K =',023.16) 


1021 FORMAT(212,214,D8.3) 

MmeoO FORA? (F1l] .7,F12.75,3X, 3018.10) 

eo) FORA tr ll. %7sF12.f+39%,D018.10) 

owe FORMAT (*'-—' .5X,*EIGENVALUES OF H FOR ALPHA = ',Fl2<7y/)} 
3004 FORMAT (4X,5D25.16) 


Poo oerOUnRha! ('-—',/5,5X-'’ THE VALUES OF K FOR ALPHA =!',FlZe/,' 
Pen «le O'R 


<ND 
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SUDROUTINE TBLI(MAX,A,.8) 


feet GENERATES THE TABLE OF BOUND-BOUND INTEGRALS 


IHPLICIT REAL*¥8 
DIMENSTON A(1),B(1) 
five INTEGRALS ARE STORED 
PO@MOING TO THE VARIOUS V 
MAX2 = 
LIM = 


IN A 
AEUIES 
MA X3MAX 
GHAX — J )*HAX2 


Meer tALIZE THE A AND B ARRAYS 
DO 1 F=1,LIM 


A(I) = .9D9 
1 @@iJ) = .O0D0 
NI4AX = MAX - 1 


Svea tTeE tHe INTEGRALS FOR MU 
DO 2 NiN=1l.NMAX 
DO 2 LL=2,MAX 


(A-H,O-Z) 


EO OLPMENS 0. AosCURRES— 


NU, LANDA 


ON AL ARR 
OF AND MU 


IFC OMWNF+LELOLT.5).0OR. (NN4+LL.~GTeMAX+2)) GO TO 2 


IN = NN + MAX¥(LL 1) 
CATE SUMCTUN, LL AA) 
AC(IN) = AA 


2 C@NT INUE 


tu Za e 


Oeomee RECURSION RELATIONS TO FIND THE INTEGRALS FOR MU = 


DO 3 NN=2-MAX 
DO 3 LL=2,MAX 


PE ( (CO LL. GT. WAX+3 ).OR-C(NNFLL.LT.G)? GO TO 3 


Th = Se MA XS (LE -— 1) 
INP = IN + MAX2 - 1 
ACIN) = ACINP) 


+ MAX2 


emaiiteee . 3). ORI NN.LT eo). ORs tNNFLECET 2.6) 70R.(LL GES 


INP = 
BC IN) = 
3 CONT TRWE 
MMAX = MAX —- 1 


IN + 
Aime )/ 6 Ud 


MAX*® (MAX + 1) 


mee 


1 


J) 





NLMAX = MAX - 2 


Meee the RECURSION RELATIONS TO FIND THE INTEGRALS FOR MU > 2 
DO 5 MM=4,MMAX 
DO 5 NN=2-eNLMAX 
DO 5 LL=2.NLMAX 
TPCCNN4LL.LT.5).OR.(NNtLL4+MM.GTOMAX+5)) GO TO 5 
ie Ni + Maveet LL - 1) + (MM — Lie tAXx2 
P= Nee a2e0MAX2 —- 1) 
INL = 2 ae aX 1) 
BWNL = IN — MAX*® (2*MAX - 1) + 1 
ACIN) = ACINN) + ACINL) - .2D1*BCINNL) 


Leer. o TeAX—-3).0R. (NNtLL.~LT.6).0R-(LL-~LE. 
F 


Bie) = (BGPNN) + BOMNL) — .2D1*ATINNL))=D 
LOFLOAT (MM+1) 


5 CONTINUE 
RETURN 
END 


2).0OR- (NN 


). 
Oars: 
iy AGU or) 


SUBRDUTINE SUN CRU,LMDA,CUT ) 


SUM EVALUATES ACNU.s.LAMDA,2) FOR GIVEN NU, LAMDA 
meePLICIT REAL*8 (A-H,O-Z) 
EXTERNAL FACT 
A = Q.0D0 
NN = NU - 1° 
DO 10 1=1,.LMDA 
IN2 = LMDA - I 
INL = NN + IN2 
eee + Ener clN 1 )/ CRAG TILING I=UO.2Drer Gl Ni to) )3 
tL = fee - 1 
meme = FACT (LLI*=CFACT(NNI — Ad 
RETURN 
END 
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De@WeLE PRECISION FUNCT TON FACT IN) 


Peet FORMS N FACTORIAL 

meerere}) REALS 8 (A-H,O-Z) 
FACT = .1D1 
tee. EO.0)) CO 10-20 
DO 10 T=1,N 
Fl = | 

Pemeeac] = FACT=FI 

20 RETURN 
END 


pity 





SUBROUTINE SPLINL(XsYsMeXINTSY INT) 


PAV EOPVALUESOT THE VGRbi NATE USING 


(a) = 


Cyc. 
~IOoO ZO 
~<—I Fe Ow 
Iormw Q. 
COLO AMO 
wast KK WU 
I VW WE 
2 ee oct 
ZIreoZz OA 
We) Oe) 
NZeO a. 
| ys ol 
MIO 
<{t NCAZ me 
od ro @ Ga peas YB Te ea wipe 
<I I NW se 
me  &— OSNY YU 
2uI27 2eusa OM 
pei oe nt OO LL 
wee YHNeidwWw 2) 
eYAe FOr tbe 
ll on Soles 2 See) 
ete, OL 7) a~aNY 


be owt ROT 
DO ee 


ORDINATE 


COR EAT RAP OLA eb) 


moe Z2iI<cf 


OY extuwxKl dite 


DHMNwZware | OO 
eee 


YN NRPS OY 


Ade de OZM SY 
<tW E Td 
OO <{ 7 <LQ. ee ee — 

OCIOHK r= ~x< 
ke Ww 


2s) ee Cee 

aI tLOOt~ MY IO FR OI 
Oogqgioruw Otazxr7a 
FAO OD (Loewe > 


YIN: 


WATHMEMAT! 


REFERENCE: 


ih) ROO Cre 


Ley SS 5 
CA lh) sh eel 


(O=27) 


aweelLicil REALS 


X(M)+YC%),C0(4,6) 


CraelL SPielC Oley. h,C) 


Perens 1 GN 
K=1 


SPLIMSICUXs YeMe XINT,YINT) 


3 “test =xX01)) 


70 K=1 


ENTEY 


Ose, 2 


COsero, / 


ret 1) 


RETURN 
2 IFC(XINT—-X(KHt1) )6+4,5 


4 YINT 


1 YINT 


eier } 


RETURN 


> ik 


=K+1 


11 satis: 3 


litetatl—K ) 


i aia 


Ge 0 7 
Osi mmeerNT—XCKI)135,12,11 


LZ Y Leal 


71 sK 


CR } 


RETURN 


Fale: 





is K-K-1 
GO TO 6 
TDW Ol, XINT 


POPeeemeret( GHOXAINT = £€18.9,32H, OUT OF RANGE FOR 
MerekPeLAT ION) 


Pema Oe J SXINT) (C1 1s K)*( X( K+1)—XI NT) ** 240 (3,K) ) 
Weerey—VINT+t CAMNTA—XACK)I*(C(2,K) F(X INTHACK) ) = 324C(4,K)) 
RETURN 
END 


SsOeroUroNne SPLICO(X,Y,4,C) 
mreemerc ity REAL*8 (A-H)s+REAL*8 (OG-Z) 
| ity ial XM) s¥(M) 601456) ,D006),P(6),E(6)¢A(6:3),B(6)- 


MM=M—] 
DO 2 K=1.MM 
D(K) =A K+1)-X0K) 
P(K)J=O(K)/6. 

2 ee ={ YC K4+1)-Y(K}) /O(K) 
DO 3 K=2,1M 

2 6G) -EUR)-—E(K-1) 
Wail 2) =~ 12-001) /D02) 
A(1.3)=D1)/D(2) 
fot 2,3)=P(2)-—P(1)*A(1,3) 
Pabst) 26. lot Pie =e tl) ae 2 
Mees S)=hWl253)/ AC2, 2) 
B(2)=B02)/A( 2,2) 
DO 4 K=3,MiM 
ieee = 2 et OS — 1) +P (CK) )SPUKRH1) -A(KH1 33) 
Pai =S 0K) =P tK-1)+8(K=1)) 
A(Ks3)=P(K)/AK,2) 

4 B(K)=B8(K)/A(K, 2) 
Q=D(M-2)/D(M-1) 
Mais 1 )=1.4+0+AtH-2, 3) 
PGs 2) =-Q-AlH, 1) =A(M-153) 
Perm). =BiM=2 )—AlH+1)=8(M-1) 


eS, 





7{M)=B(M)/A(M, 2) 


MN=M—2 
DO 6 T=1,MN 
K=M-] 


Z(K)=B(K I-A Ks, 3) ¥Z(K+1) 
Merah) <2 (2 )-Al(1 +3) *Z(3) 
DO 7 K=1,MM 

O=1./(6.*D{K)) 

C(1¢«KI=Z(K )*Q 

Cl2sK)=Z(K+1)*Q 

Os SC K) /DCK) -—Z(K) =P (KK) 

Oe MK I=HY(K+1)/0(K)-Z0K +1) *PUK) 
RETURN 

ERD 
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SHB me) INE ELEMGAH, SS+ Ay By MAX) 
PicdI(d}> FOR HAND S, RESPECTIVELY 
Prec l@dlT REACeS (A-H,0-Z) 

INTEGER PWR 

COMMON/TRM/ALFA, ETA 

COMMON/ATNO/Z 

COMMONS INDEX/NIYNJIsLIsLJeM1,MIS, 1S 

DEMPENS ION Ati) .8(1) pNU(52) pAU( 52),LMDA(52) ,COEF(52) 
Niwe=-eN!] + Ne 

ie = LI +L 

Nerd =ONT + °C 

Caw = LI + NJ 

i) re a ee aS, 

FN1 OFLGAN I + 1) 

FLL ODELOATI{. J + 1) 

FN = DFLOAT(NJ) 

ke OFCOM TCT) 

FM OFLCAT(MJ ) 

MAX2 = MAX*MAX 


ES OF UNS, LAMWDA, MU ANDO TH 
D2 | Ghee OF HIl,JIo AND THE 4 


NU(C 1) = NIJ + 2 
NEE 3S) = NU (de) 
NU( 4) = NU(1) 
NU{ 6) = NUCL) 
NU( 8) = NU(1) 
NUC 9) = NUCL) 
NU(37) = NU(1) 


m 
O 
ak) 


LMDA(10) = NUC 1) 
LMDA( 11) = NUC 1) 
LMDA(14) = NUC 1) 
LMDA(15) = NUC 1) 
LMDA(1E) = NUC 1) 
LMDA(18) = NUC 1) 
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LMDA(42) 
NU( 2) = 
LMDA€12) 
NU(39) = 
LMDA(44) 
NUC 5) = 
NU(38) = 
NU(40) = 
LMDA(13) 
LMDA(41) 
LMDA( 43) 
NU( 7) = 
LMDA(17) 
NUC 10) = 
NU(12) = 
NU(13) = 
NU(15) = 
NUC17) = 
NUC18) = 
NU(41) = 
LMDA( 1) 
LMDA( 2) 
LMDA( 5) 
WDA ( 6) 
LMDA( 7) 
LMDA( 9) 
LMDAL 3 8) 
NU(11) = 
LIDA( 3) 
NU(42) = 
LMDA(40) 
NU(14) = 
LMDAC 4) 
NU(42) = 
NU(44) = 
LMDA(37) 


= NU( 1) 
Nid + 1 
NU( 2) 
NU(2) 
NU( 2) 
WJ + 3 
NU (5) 
NU(5) 

= NU( 5) 
= NU(5) 
= NU( 5) 
NIJ 

= NU( 7) 
rom one 
NUC19) 
NU( 10) 
NU(C1)) 
NU( 10) 
NU( 10) 
NU(19) 

= NU( 10) 
= NU(1)0) 
= NU( 10) 
= NU(10) 
= NU(10) 
= NU(10) 
= NU(1)) 
oe eas co | 
= NUC11) 
NU (il) 
=—awWile 1 ) 
eile: 3 
= NU(14) 
NUC14) 
NU(14) 

= NJU(14) 


Ns 





LMDAT39) 
NU(16) = 
LMDA( 8) 
WEIS ) 
NU(21) = 
NU(22) = 
NU(24) = 
NU( 26) = 
NU(27) = 
NU(45) = 
LMDA(28) 
LMDA(29 ) 
LMDA( 32) 
Emer 3 3 ) 
LMDA( 34) 
LMDA(36) 
LMDA(50) 
NU(20) = 
NU(47) = 
GeDA( 30) 
LMDA(52) 
NU(23) = 
NU(46) = 
NU(48) = 
LMDA(31) 
LMDA(49) 
LMDA(51) 
NU(25) = 
WMDA(35) 
NU(C28) = 
NU(3Q0) = 
NU(31) = 
NU(23) = 
NU(35) = 
NU(36) = 
NU(49) = 


it 


= NU(14) 
oS 

= NU( 16) 
Sibu 2 
NU( 19) 
RUE ne: )) 
NOES ) 
NU(19) 
NU(1)) 
NU(19) 

= NU(19) 
= NU(19) 
= NU(19) 
= NU(19) 
= NU(19) 
= NU(19) 
= NU(19) 
els +. ] 
NU: ) 

= NU(20) 
= NU( 20) 
wc TS +. 3 
NU( 23) 
NU (23) 

= NU(23) 
= NU( 23) 
= NU(23) 
NLIJ 

= NU(25) 
Pivot 2 
PolrZ Ss ) 
NU(28 ) 
NU(28) 
NU(23) 
NU( 28) 
NU( 23) 


L253 





_ 
ae 


fe eFF0e Fee 
it 


a 


a 
— 











LMDA(19) 
LMDA(2Q) 
LMDA(23) 
LMDA(24) 
LMDA(25) 
PPD ANE 2 7") 
LMDA(46) 
NU(29) = 
Newol) = 
LMDA(21) 
LMDA(48 ) 
NU(2Z2) = 
Smet 2 2 ) 
NU(50) = 
NU(52) = 
LMDA(45 ) 
LMDA(47) 
NU(34) = 
LHDA(26)} 
MU 1) = 
MUC 2) = 
MU( 3) = 
MU( 7) = 
HU 8) = 
MUC 10) = 
MU(11) = 
rovrc) = 
MU(16) = 
MUC17) = 
MUC19) = 
MU(20) = 
MU(21) = 
MU(25) = 
MU(26) = 
iU( 28) = 
MU(29) = 


NUT 28 ) 
NU( 28) 
NU( 28) 
NU( 28) 
NU( 28) 
NU( 28 ) 
= NU(28) 
ears + ) 
NU(29 ) 

= NU(29) 
= NU( 29) 
ee + 3 
= NU( 32) 
NU( 32) 
NU(32) 

= NU(32) 
= NU( 32) 
L Nine 

= NU( 34) 
mis + 2 
MU(1) 

MU 1) 
MU(1) 
MU(1) 
MU( 1) 
MUC1) 
MUL) 
MU(1) 
MU( 1) 
MUCL) 
MU(1) 
MU(1) 
MU(1) 
MU(1) 
MU( 1) 
MUC1) 
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MU(30) 
MU(34) 
povs 5 ) 
MUC 4) 
MUC 5) 
MUC 6) 
MU(13) 
MU(14) 
MU( 15) 
Pree? 2 ) 
MU(23) 
MU(24) 
mo 3) ) 
MU(32) 
MU (33) 
MU(37) 
MU(38) 
MU (39) 
MU(40) 
MU(41) 
MU(42) 
MU(43) 
MU(44) 
MU(45) 
MU(4C) 
MU(47) 
MU(48) 
MU(49) 
MU( 50) 
MU(51) 
MU(52) 
wit 9) 
MUC18) 
ren 2 i) 
MU (36) 
COEF { 


1) 


Midi) 


MU(1) 
MU(1) 
MIJ 

MU (4) 
MU(4) 
MU(4) 
MU(4) 
MU(4) 
MU (4) 
MU(4) 
MU(4) 
MU(4) 
MU 4) 
MU (4) 
MUC4) 
MU(4) 
MU C4) 
MU(4) 
MU (4) 
MU 4) 
MU(4) 
MU(4) 
MU( 4) 
MU (4) 
MU(4) 
MU(4) 
MU(4) 
MU( 4) 
MU(4) 
MU C4) 
Mel Jet 
MUD) 
MU(9) 
MU) 


if 


= —0. 25D0=ALFAALFA 
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aid 


COERO 10 
Cwer (19 .) 
COSPY 2 8°) 
COBF( 2) 
C@ireret rz) 
COEr Tay) 
COEF (Z9") 
CGIEF({.-3) 
COEF (1.13 
COmry 2 0") 
COEF (3 98 
Gere { 4) 
Cees) 
CGEr (13) 
CGEF(14) 
C@eEl 22) 
COEF (23) 
C@Er( 31) 
Oger 132} 
Geerror} 
COEF (3S) 
COEE(41 ) 
COeEA42 3 
COS (45 } 
COEF (46) 
COEF (49) 
Ger (50) 
Geer (G.) 
COEF (iS) 
COEF( 24) 
beer( 43} 
C@Er{ ££) 
COEF(17) 
COEF (26) 
COEF(34) 
GCOEFAn«8:) 


COEREh)) 
COEF()) 
COE Rta) 


J-5 DO AL FA* FNI 


CUSmHac? 
CWE C2) 
COEREZ ) 


Ose RO ALFA=FLI] —-Z 


COREG) 
COE RRS) 
COERT3) 


0.5D00*ALFA*FM 


COEF(4) 
COEFFI) 
COEF( 4) 
CHEE (4) 
COEF U4) 
COEF () 
COEEL4) 
=wer (4) 
=eoer (4) 
=60EF (4) 
Se@rr | 4} 
-~CUEF (4) 
=COERC4) 
—-C@EF (4) 
—-COEF(4) 
dle irl oe whl 
CUEFTG) 
COEF (Cc) 
GOEF (ey 


—0, .00* FN AEN] 


CORP 7) 
COERC?) 
Gwe (7) 


=WeoO0s Fletr LL 


5 et | ol 





O@er (TG = COEF(8) 
C@emr25) = CO=F( 8) 
Geer (359) = COEF(8) 
G@eE{ S).= 0-101 
CHER( 16), = -COEFI9) 
Vbena,}, = COEF(9) 
Ceer{ 36a = COz=F(9) 
oeer(3o)> = FRFFN 
COEF(44) = COEF( 39) 
CGmat43) = COEF(39) 
Cie 51) = COEF(39) 
OBEF(40) = FM*FFL 
COEF(43) = COEF(40) 
COEF(47) = CCEF(40) 
e@erto2) = CUEF(4)) 
MAX2 = MAX*MAX 

HH = 0.0D0 


Herts .20.0) 


Gee 


GO Ste sta UND Wars eles scar 


POR WAIPLET 
med. J) MUST BE CHARGED 
DO 1 I=19.36 
feoriet li) = =COERET ) 
DO 2 1=45,52 
peeetemi! ) —-—CHEr( 1) 
Poe aVSiNeseinbex ESMRESA 
INTEGRALS 
10 DO 6 K=1,.36 
Pree NURS + 1 + MA 
PWR = NUC(K) + LMDA( 
6 HH = HH EO CUE 
DOE POK=3 7452 
IN = NUC(K) + 1 + MA 
PWR = NU(K) + LMDA ( 
7? HH = HH 
eae yy DE tJ) 


ND 


ONDING 70 THE STORAGE 


X*LMDA(K) + MUCK) *MAX2 
K) + MUCK) 
F(K)*AC IN) ALFA** PWR 
X*LMDA(K) + MU(K)*MAX2 
K) + MUTK) 


WR 


Mie. S6heePROCE sores BE 


SCATIERMWG THE SIGNS OF THE LAST 26 TERMS IN 


MU FOUND ABOVE MUST BE CONVERTED 
NODE Oe wilitctt 


REPEATED 





1) 


Ll 


O18. J. 0D0 

DO 9 K=1,10-9 

Pee NUK 5+ 1 EoMAXSLMDA{K) + MUCK) *MAX2 
pyos = NUCK) + LHDA(K? + NUK) 

wom= So + Ati Air A*= PIR 

Weer is.60.0) GO TO ll 

De We K=19,28,9 

IN WUTKY + 1 + MAX=LMDACK) + MU(K)*=MAX2 
PWR = NU(K) + LMDA(K) + MUCK) 

ae SoS ACING7 ALF Ax=PWR 

GO TO 12 

DES S K=1:9..28,9 

fee INC) + ) UE MAXSLMDA(K) © MUCK) MAX 2 
Pre = WICK) + LMHDA(K) + MUCK) 

wor oo t ATINI/SALFAs* PWR 

RETURN 

END 


It 
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SwereOorINE SCL{(S+NyF) 


OF The SEARCH Sd (ABSOLUTE 


> 
Or Ge 


el TCA L * 

DIMENSION $1) 

TOP = DABS(S(1)) 

DO 1090 T=1,N 

T = CABS(S(1)) 

tee tieGi=zh@P) TOP = 7 
100 CONTINUE 

Gee= DLOGLO(TIP) 

IPWR = IDINT(FF) 

pen eee LP) 

RETURN 

ENO 


SUBROUTINE SM2°YCA,C WN) 


URNING THE SCALED 
AL MATRIX 


<4 


mw LICIT REAL= 
DIWENSTON Al) 
Be & lotta 

1 ALT) = ATIC 
RETURN 
EMD 
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CA iis Oa 


SUBROUTINE INVERT(TSsH+S+C.N,ISIZE) 


CALL DSINV(TS,CyNyEPS,TER) 


iemlER.EQ.—1 ) 


MeENSTON 1$(1).H01)+801),CCO1) 
WSelG 


WaPLICIY REAL*8 
KIT 


REAL EPS 


GO TO 16 


2aN 


DU 3 K 


eZ 


= K¥(K - 
Mor2 1=1.KK 


Ki4 


DO 1 J=1-K 


ee eg | 


d= 


— 
— 


IJ 


CC + CCl J = T SG) 


ee 


0 


— 
me 


IS1G 


=e! Oe 


2 we) 
3. Ceny MAUVE 


4 
B 


Sie Gi 


MAIN PROGRAM CAN 


DY GOOD ENOUGH 
GO TO 15 


E 


<I 


IF(ISIG.EQ.1) 


14 KIT 


KIT + 1] 


Me TERM IS FORMED BY MULTIPLYING U**(-1) AND THE 


DC 6 K=2,N 


0) 








ee= K — ol 

ieee Kol — 1)/2 

DO 5 IT=1+KK 

IN = Ne(I — 1) 

IK = K + IN 

Pie= >] + 1 

CC = 0.000 

CO 4 J=II-K 

Igie= J + aN 

wee = WH + J 
face] CC em rel) *TS (JK) 
> Stile) = CC 
6 CONTINUE 


en MeERSE 1S FURMED BY ADOING THE CORRECTION TeR™ 


DO 8 K=2,N 

Mae Ko =.) 

he = KeeK -— 1)/2 

Dory I=lsKK 

IK = KM + J] 

KI . eee 1) 

Wwamtk? = TS01<) + S(KI) 
S “CON TANUE 

ISIG = 1 


THE NEW TRIAL INV 
> 108€(-20). THI 


DO 11 K=2-N 


ne ae 

DO 10 J=1.KK 
IN = Neal. = 3a 
a en 
CC = 0.0D0 

If =] +1 


PREV iwGl chim) Somme | 0 
DO 9 J=I11-+KK 
IJ = J + IN 


Boa. 





1) 


Cesta) = ta St) 


oC. = 


= 0 


Won G 


Tf(DABS(CC).GI.0.1D-20) 


hee strh) 


Ce 


11 CONTINUE 


CON 105 


Mp AlesiG 6C%.1 } 


Hew TRA TN= 


2aN 


oe is Kk 


=1,KK 


per r2 | 


1) 


ee ey lee 


[Kes 


LSeGOMT. have 


GO TO 14 
SLT E (64,1001) 


oe. 


RETURN 


WRITE(6,1000) 


sep 
1000 FORMAT (10X.'STOP. 


16 


AAO SeGaninO) BE FACTORED.*) 


=a 


Ee 


BOOlL Feat ('="',1GX,"* INVERSE FOUND AFTERSs12;'ITERATIONS.*) 


END 


Svor.0UulMee DSINV(AsC .N EPS, IER) 


bid tS) 
VY) ere a 

ce —Y (Bs Nz 

we) ieee eas ZZ HOE KH 
<> St a =) 

— qa oy > Te Oo 
eau kb Ut me OW WwW 

ke C0 OTK oO Oo eu! 

Dal T— On & Th- WN 
re Ge Oe Oe SG ioe «I 

ea = 8 ae Be Ee 

WU Miva. Ee EO. 

k- > Ire = oO ee (il 

— J pee | ome em ee. aS 

Wee ee Dal 9 A, 9 Ge oe ced © 
—O (Ory <I | Se aes 
ul. grates) (Cle es 2e20 

Lj UJ JireD KY NN Ox ce 
ae Mm Ow EOS te te Re” tL 
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EES: 
IER 


Because OF WRONG SNe 
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Pon eee eon. DECAUSE Sone RAD 
PEANO VS gINON=POSDLLIVE {MATRIX A 
Poo Poot Vester Ui es POS 
Sie roves oe WS > Geess| GNirfi— 
CANCE) 

TER=K WARNING WHICH INDICATES LOSS OF 
Stonir fCANCE. ~ IHE RAGTCAND 
BeICIED Aly *FACTORIZATION STEP Ket 
toe DE ePUS i 1 VesbuT NO LONGER 
Pees AKL iy itl) | 

DOUBLE PRECISION AsDIN »WORK,C 


DIMENSTON A(1),C(1) 


mo THE WPPER TRIANGULAR FACTOR OF A 
GAL DWF SD(A.N,sEPS,IER) 
MeClER) Sy1lat 


INVERT UPPER TRIANGJLAR MATRIX T 
Pmereee IfiVERSTION LOOP 
1 IPIV=Na (N41 )72 
LND=1PAiV 


Seeree «6 CY DUPLICATING IT IN C 
D@e/ T= IPIV 
7 Gl) = AT) 


ieee i Ze -PNVERSTON-LOOP 

DO 6 I=1.N 
DiN=%.CD0/ACIPIV) 
A(IPIV)=DIN 
MIN=N 
KEND=I-1 
LANF=N-KEND 
BROCKENG ) 54542 

2 J=IND 
DO 4 K=1,KEND 


PRITEALIZE RCW-LOOP 
WORK=0 .DO 
Pb =M Tilt 1 
LHOR=IPIV 
LVER=J 


START INNER LOOP 


Le, 





DO 3 L=LANF,MIN 
EYVER=LVER+1 
LHOR=LHOR+L 
3 WORK=WCRK+A(LVER) *XA(LHOR) 
A( J) =-WORK=DIN 
J=J-MIN 
IPIV=IPIV-MIN 
IND=IND-1 
RETURN 
END 


0 & Wn 


SWBROUTINE OMFSO(A,N,EPS,IER) 


OLMENSTON: Al(1) 
OS aee, PRECISION DPIV; DSUMRA 


best GH WheeNG TMPUT PARAMETER WN 
Memin—l) I2e1 yl 
1 IER=0 


Met TALLTZE DIAGONAL LOOP 
KPI V=9 
DO 11 K=1,N 
KPI V=KPIV+K° 
IND=KP 1V 
LEND=K-1 


CALCULATE TOLERANCE 
TOL=ABS(EPS*SNGL(A(KPIV))) 


Sue. FACTORIZATION LOOGP OVER K-TH ROW 
Dell =Kh yh 
DSUM=0.D0 
eNO 2 6 2 
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START INNER LOOP 
2 DO 3 L=1.LEND 
LANF=KPIV-L 
LINO=1 ND-L 
3 DSUM=DSUMFA0LANF) *A(LIND) 


TRANSFORM ELEMENT ACINO) 
4 DSUM=A(CIND)-DSUM 
ie i—K) 10s 5510 


eae Oe WEGATMME Plaey ELEMENT AND FOR LOSS OF SIGNIFI- 
See coNGLADSOm)-10L) 6:649 
6 IF (DSUM) |p aes | 
TVIPCIER) 85,9 
Gr PERK =K- 1 


COmPwTE PIVOT ELEMENT 
9 DPIV=DSORT(DSIM) 
ACKPIV)=DPIV 
ODP IV=1 .DO/ DP IV 
Cert0 if 


WAPCULATE Terms IN ROW 
10 ACIND) =OSUM*DP TV 
11 IND=INO+F] 
RETURN 
pe] LER=~-1 
RETURN 
END 
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SUBROUT WWE TMPRO(HU,CyNe F) 
tee FORMS THE MATRIX PRODUCT U* 
eee eile > 11 BY THE SAME FACTOR 

ee ol REAM SC A-HyO-Z) 

DPMENSION HCl) ,UMl) @¢ 1) 

Ea = NetN + ol) 2 

DO 3 L=l.N 

DO 3 K=L,N 

LK = Ketihee li72 + 

UHU = 0.000 

DO 2 T=lel 

HU = 0.000 

Ver = Eek = Dees: vl 

DO 1 J=1.K 

i= Jar) — liege | 

Mae GTleJ) 1) = I¥(I - 1)/2 + J 

Sikte= (Kanth =) ieee oJ 

LeAU = HU + HC IJ) USK) 

UHU = UHU + UCILI*HU 

3 CULK) »= a0nRv 
DO 4 IT=1.tLI 
ye) = °C Clear 
RETURI 
END 
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SUBROUTINE DJACVT (CA Ns, NOYES,EI VU,EIVRsNDIM) 
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Wk meer uO 
D DW Wort Te 
ee a ee | AD 

wo raze weke 2 
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ZOwwWY fF 


eM WE) 


Or eCwrvOU>S>N 2ZOOWrt WeEZUWOOWY OWZWOUW 
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Re SIZE RY ENO OM DE SITE 


= eS ae MU ae 
Sr“ ZOU OT O ee 
CL SL Ter WY 


Mk © OFrRWOOUW 
ULI eo (ee > 
k= OO CM wae ree) 
See ee LEEY. 7 Le. 


Eley: 
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~ <[ Y) 
2, eyoea) i 
e f= SLL =a 
saa bw Coe Oe 
> kl 
me ee YL 
b—b- LW co 
<I Wiis 

ee C tase ae | 
DED) ae) =) 
Oat <n as “> 
JW kw © 
<{ > Ww oO 
CG) we 

lw ZR eu = 
Or ee 
2<%I WN> DE 
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me | 3 Gee Fe 
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Oc Wid 
uJ = OQ ee N 
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Sek) Core <t 
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Sood © 
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OY OnHAtInM Ss 
uJ @WooewoO ae a 
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a I ase J ED | es 1 
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METH GD 


JOWRNAL OF AGM. 1, 


eee c8.5 | 


Si!, 


~U 


Or 


meee ReAL “Syaam 
1359 


WILK IW 


S 
VON NEURAN, Jeo AN 


Ret ere CE 


PROG Ee 


ON, Je Hes 

ORD UNIVER 
PeePLICIT REALS (A-H,0-Z) 
DIMENSION ACNOIM,NDI 


lA =—1)20,23,.21 


20 PRINT 22,N 


4 
t 


M) sclLVUCNDIM) »EIVRINDIM,NDIM) 


RETURN 


ei lo? le 


RETURN 
23 PRINT 24; 


Al; 1) 


= eel aad} 


MATRIX A 


. 


24 FORMAT (24H IN DJACVT 


RETURN 


a7 





2) WEG Sh6.9)1+ly3 
3 PRINT 2,N 


é FUORMMAIC' wW= ',13 
LPOePCACLING ROJTIN 


RETURN 
1 IF(NOYES)99,192,99 
99 CONTINUE 
DO 101 J=1,N 
DO 100 I=1,N 
100 EIVR(I1,J)=0.0 
oie ¢ J,J)=1.0 
102 ATOP=O0. 
DO 112 J=1lsN 
DO 111 I=leJ 
i SUMO EEW IRIE Vere auiee scle 
99 PRINT 106,NsN 
TOG FORMAT (14H IN DJACVT (Al,13+1H,+13+3H))9) 
PRINT 108.1 .J, J; 
; 


MOO ORMAT ment ta) oat yg 'gi3y') AND Atty 13¢',"', 
PRRNOUAL. THEY WERE REPLACED BY THEIR MEAN 


Rete dd= .<Se0ACI,JI+A(J,1)) 
A(J+T)=A( I,J) 
103 CONTINUE 
PP CALGP [Oe estA Cle Jd) 7) 104s lives lit 
104 ATOP=DASBS(A(T,J)) 
111 CONTINUE 
Liz BIVU(JI=ATI.J) 
eaters yl I9.113 
L109 ©P Rail 
PGR SORMwee26H IN DJACVT, MATRIX A = 0 } 
RETURN : 
L113 AMGF=ADFLOAT(NS(N-1))*.55 
D=J.0 
DUP Tre Jd=ZyN 
DO 114 [f1=27JJ 
S=ACTI-1,JJ)/ATOP 
114 D=S*S$+D 


Lame Wes: Bit toto oo. RETURN 
Sour Reh 
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117 


Ls) 


116 


120 


EZ) 


127 


1 2e 
23 


Dor El 0 =Car 0 

Peesh =DSORT (D7ZAVGF)*ATOP 
IFLAG=0 

DOGO. JCOL=2ak 

Jeol 1=JCOL=1 

D@r 130» 1R@W=1,9C€0L1 
AIJ=A(TROW,JCIL) 

Peveno stay J) -TARSH) 130.150;117 
ALT=A(IROW, IROW) 
AIN=A(JICOLwJC 3b.) 

S=AJJ-AIl 

erenees (Ai J)—-1.0D-09*DABS(S))11304+130,118 
IFLAG=1 

Drees D—- 1O0"DABS (ATJ)-—-DABS(S111164119,119 
S=. 707106/7811865475 

C=$ 

Go, 1B»120 

T=AIJI/S 

$=30,. 25/DSQRT (0. 25"T*T } 
C=aDSORT (0. 5+S} 

SZ sola C 

OOmra T=). haan 

T=A(1,TROW) 

U=A(1,JCCL ) 

ACI ,IROW)=C*T-S*U 

AM)» JCC’) =S= Tec *U 

I2= IROW+ 2 

Pel 2—= JCC i Z2islZisl2s 
CONTINUE 

DO 122 T=l2,J-Q0L 

l2aGiaie COL) 

U=A( IRCW,I-1) 

SO — 15 SOOM) =S*U+C=7T 
ACTROW, 1-1 }=C¥JU-S*T 
Adm@lCOLs JCOL) =S*Al J+ C=AIJ 
ACIRGW, TROW)=C¥A(TROW, TROW)-S*(C*¥AILJ-S*AJJ) 


Loo 





DO 124 J=JCOL,N 
T=A(ITROW.J) 
U=A(JCOL,J) 
A(TROW,J)=C#T-S*U 
124 A(JCCL»J)=S*T#C#U 
TF(NOYES)131,126,131 
131 CONTINUE 
pe 125° PSP N 
T=EIVR(I,IROW) 
EIVR(I,ITROW)=C*¥T-EIVR(I,JCOL)*S 
125 EIVR(I,JCOL) =S*T+EIVR(I,JCOL)*C 
126 CONTINUE 
S=AIJ/ATOP 
D=0-S*S 
IF(0-OSTCP)1260,129,129 
paso = 0). 
OO.128 JaeeeN 
em 26 (lea 
SSM 1I-1,JJ5) 78T 0P 
128 D=S*S+D 
DSTOP=(1.0-06)*D 
129 THRSH=DSORT(D/AVGF) *ATOP 
130 CONTINUE 
IF(IFLAGI115,134,115 
acy T= Aes) ) 
AGS FIV) 
EIVU(1)=T 
DO 132 J=2,.N 
WARIS) 
We, J) = 2 IV) 
El VEuGak= 7 
Oe a2 5-23 
Pe 2 esti A (el 1) 
133 RETURN 
END 


140 
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SUBROUTINE MPRD(A,B,C,N) 
feeeo eens JHE MATRIX -PRODUCT A*Bs JHE RESULT IS RETURNED 
Meeae C [S$ USED FOR TEMPORARY STORAGE. 
ieee cir REACTS CA—H,0=Z) 
DPE NSTON AGL) »B'h)sCO1) 
DO 2 1=1,N 
DO 2 K=1,N 
wm = N#K-1) + I 
CC = 0.0D0 
DO 1 J=1,N 
Pe = Wet J—-1) + | 
JK = N¥(K=-1) + J 
gee eG = CC te Al ld) eB CIK) 
Peek) = <CC 
Lift = NN 
DO 3 1=1,LIM 
Sree le = (CCl) 
RETURN 
END 


141 





SUBROUTIWE INI 256( FCT, F) 


WMELICIT REAESS (A—H,0=Z) 


COMMON/L256/X (128) ,A(128) 


COMMON /CD/O0PC(15),;D0MC(15) 


0. O0D0 


DO 100 


Lye 6 


— 
— 


=1,15 


DO lov J 


= DPC(J) 


OP 


DMC(J) 


DM 


Site ben VAL RES 


CORE ORT 10. LHe 
KE QUADRATURE 


OO 


XX*DM 


LLL. oe CTP A Pei) 


Fo + ZZZ*AA 


MOP 227 


100 F 


RETURN 
END 


DOUBLE PRECISION FURECTION CGARG(Z) 


AND AT THE DESIGNATED POINT FOR 


at 


CA-H,0=—Z) 


Pee Cl Tt eae 


Va 





COMman/CONST/PH1,PH2,FR1,FR2,SCL,PHOP] ,FCT2,FCT1,A 

ARG = AX®DLOG(Z! 

Oe Ge i aes) VEX (—- 2 Cl 2 OS IN CARG) + FCTL*OCOS(AR 
RETURN 

END 
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SwerOUTINE TBEZCIMAX, MAXI; AS,+AC; BS,BC) 


TBL2 TABULATES THE REQUIRED BCUND-FREE INTEGRALS 
IMPLICIT REAL®8 (A-H,O-Z) 
COMMON/CONST/PH1,PH2;FR1yFR2,ASCLsPI¢CArCTsAAA 
COMMON/WAVP/WAV 
COMMON/TRM/ALFA, ETA 
COMMON/ ATNO/Z 
COMMON /RHO/N 
DIMENSION AC(1)+AS(1),8C(1),BSC1) 

DIMENSION A1(144) ,A2(144) 

DIMENSION $1(13),C2(13),SS(11),CC(11) 
DIMENSION X0(500) 

ALFA2 = ALFA + Z 


a ie Pn eeeenC Yon PHASe FACTORS IN THE SINUSOIDAL 


FRILL = WAV/ALFA2 

PH] = mea DLEGt.2DL=FRI) + ETA 
FR2 = .~2D1*WAV/ALFA 

PH2 = wy Pee Cl a. ZDIerR 2) + EDA 


ASCL = OQ.5DO*ALFA/ALFA2 
Mei Tete. 200) Phir Rly Phl2yFhZ sy ASCL 


IF Z=1 THE NUMERICAL INTEGRATIONS NEED NOT BE DONE 
IFCAAA.EQ.0.000) GO TO 104 


me D The ZERGeES OR THE INTEGRAND FOR 12 
CABLE ERO Cid, XO; 10) 


EVALUATE THe I] SWAT EGRALS 
DO 100 I1=1.IM4AX 
Ne= bee 
CAE Sea yt LO XO, 275 ASCE) 
190 S1(I) = ZZ 


GND TREeZer ees sO tte TATEGRAND -FOR. 1S 
ALL ZERO(2,L, X0,1Q) 


144 


ee 





SveeontlE THE IS INTEGRALS 
DO 101 I=1,.MAXI 
Mee 1 - 2 
CALL CSINT(2,L,1Q;X0,22;0.101) 
Pete esol) = 22 


mre Ime ZEROES OF THE INTEGRAND FOR 13 
CALL ZERO(3,4L+ X0,1Q) 


BeecUATE THE 13 INTEGRALS 
DO 102 J=1l;IMAX 
N = I - 2 
CALL CSINT (3+sL+1Q;X0s 22, ASCL) 
mee C2tr) = ZZ 


om OD IME ZEROES OF THE INTEGRAND FOR IC 
CALL ZERUI4,b;X%I+19) 


EVALUATE THE IC INTEGRALS 
DO 103 I=1sMAXI 
we= I - 2 
CALL CSINT(4,L,.10-.X0,22;0.1D1) 
roo CCwe = 22 
GO TO 105 


Pee 2=1 EVABUATE Whe RECUIRED INTEGRALS EXACILY 
LOS S@ALL HeaaisiGSiteG@z,SS;CC, I MAMGMAXI) 
105 WRITE(6,2001) 
MRITE(G;, 20007 MIR SI CI), C20J)sJ=1, IMAX) 
WRITE(6;2002) 
Wie tevowuzogo)(JeSS (0J),CC(J),J=1, Mar!) 
MAXII = MAXI + 1 


Diet Petr Ze THEm@eee AC, BS AND BC ARRAYS 
DO 199 J=1.1452 
AS(J) ~ODO 
BSUJ7 = 200 
AC(J) = .0D0O 
199" BEV), = SObo 
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LPT IAGEZE THE Al ANO.A2 ARRAYS 

DO 200 J=1,144 
Al(J) = .0D0 

200 A2{J) = .O0D0 
NNMAX = MAXI 
LLMAX = MAXI 
NN NNMAX + 1 
ee = 2X LX 1 
AP = 0.5D0*AL-A 


filefeAP = sQwebD1/ AP 

ALFA4 = 0.201/(0.2D1*Z + ALFA) 
USING THE IS AND IC INTEGRALS FOUND ABOVE EVALUATE ASCNU, 
LAMDA.2) AND AC(NU,LAMDA,2) 


DO 201 NN=AL+NNMAX 
SSN = SSCNN) 
CCN = CC(NN) 
AiRe = (eRewALFAP 
F = 0.101 
ALFA4P = 0.101 
DO 201 Lk=1,LLMAX 
ALFA4P = ALFAGP=ALFA4 
PeaeirOr LOA TT bly) 
leet (RELL ALT 24) - OR .(WNELLAGTSIMAX)) GO TO 201 
INA = NN + NAM¥(LL + LLM) 
TAA = FeALFASP#AP 
ASCINA) =TAA*SSN 
AC(INA) =TAA®CCN 
201 CONTINUE 
NIMAX = MAXII 
LLMAX = MAXTI 
IMAXT = IMAX + 1 
NO ThE | IANID 13 INTEGRALS FOUND ABOVE EVALUATE A1(NU, 
Av2) AND A2(NU,LAMDA, 2) 
DO 202 NN=1,N1LMAX 
DO 202 LL=1,L1 MAX 
PEt is Gl chix! OR SONNE SLs 5) GO TG 202 
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p= 
OWN 


IN12Z2 = NN + NIMAX¥* (LL —- 1) 
CALL AL2SUMN(NN, LL AASS1) 


Al{IN12) = AAG 
CALL A1L2SUM(NN.LLsAA,C2) 
A2(IN12) = AA 


202 CONTINUE 
NNMAX = MAXTI 
we) = Cran 1 


NLM = NNMAX*LLMAX - 1 

NNM = NNMAX - 1 

er = OINNEX= CO LMAX + 1) - 2 
ING THE RECURSION RELATIONS EVALUATE AS(NU,LAMDA,1), 
OU. VAMOAW1), BSOWU, GANDA, 1) AND BC(NU,LCAMDA,;1) 


DO 293 NN=H2Z],NNMAX 

EGeZOS LL=1,EEmAx 

DEAR Neer .GT oieaexi) OR. (NNtLL LIS) ) GO TO 203 
INL = NN + NNMAX*(LL + LLA1) 

INAL = INL + NLM 


INLA = LL + NIMAX#(NN - 1) 

INNA = LL + 1 + NIMAX*(NN - 2) 

AS(INI) = ASCINAL) + ALCINLA) -  ALCINNA) 
ACCLINL) = ACCINAL) + A2C(INLA) ~ A2CINNA) 


Peete Ll. eee at ELL ee) OR ste ESL a6). OR. (NN.GT. 
Died? .OR. (LL.G?.LLti)) GO TO 203 


Piibl = INL + LNM 

Web = Clo — oi +> NI MAX= NN 

WANG = LL ™+ 2 + NIMAXS (NN — 3} 

Doi Ne =) CS Noe ieee — SAL CIINNE)) /.2 01 

Be Cur = SAC Ciel) SP AZt Ieee = AZCINNG) )7. 301 
203 CONT INUE 

MMMAX = MAXI 

Waxls = TMAX + 3 

ir = TRAX — 2 

NN1 = 2*(NNMAX*=LLMAX - 1) 

NN2 = 2*NNMAX*LLHM1 

Miti3 = NNMAX=LLM -— 1 
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USING 
AG (U, 
MU 2 


204 


2000 
2091 
2002 
3000 


THE RECURSION RELATIONS EVALUAT 


LAMDA,sMU), BSC NU,LAMDA,MU) NO BC INU. LANDA cHU) FOR 
DO 204 Mi=4,M4MMAX 

Bes) = WF LOATCUMMe — 3)70FLOATINM + 7) 

DO 204 NN=2,NNMAX 

OO" 204 LL=1,LCUMAX 

IF (CNN+LL+MM.GT .MAXI3) OR. (NNtLL.LT.4)) GO TO 204 
TN = NW + NNMAX=(( LL - 1) + LLMAX*(MM - 1)) 


INN = IN — NN 


INL = IN -— NN2 

Pe = IN — NN3 

ASTIN) = ASCINN} + ASCINL) — BSLINNL)*.2D01 

pewrN)? = ACCINGP + ACTINEY — BCCINNL)*.201 
a ee ORS OR LANG Ls 
poe) = ( “BoGRNN) + BSCtNtLIi~ = ASCINNL)*.201)*FM31 
Bi ee (Boe SCUINLE) == “AGUCINNE).2D17=FM31 
CONT INUE 

RETURN 

PORT AT (4 OX eee 20). 16) 

FORMAT ISS S79 34X +! S1's24X,_' Coal) 

FORMVATUS// 7534 X- ' SS'y24X_' CO / J) 

FORMATE "+ LOX PHI = 'sh2e. 16; 7”, 1LOX,'FRI = ',023. 
lO ae Oo eee ee i. — Seer De oe OU gong TR2Z = ',D23.16 
Cae J yl eel =) 022.1647 5 "=") 


END 
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SMORQUTINE ZER OCICS +J,X0,1Q) 
mrOseaAtt THE ZEROES OF SINCAX +. B= 
eRe Be Lat X) + C} BETWEEN X=XSM AND 

Meyperctiy REAL&S (A-h,0—Z) 

COMMON/CCNST/PH1,PH2,FR1sFR2,ASCL,PI CACTI 1B 

COMMON/INT/XSMSFFeBOsXLM 

DIMENS TON X01) 

X = XSM 

Oe (0 (2ls22,21s,22)31Cs 
21 A = FRI 

CC = PH1 
cer yO 1 
le = FR 
ee = PH2 
1 ARG = A*X + B*¥DLOG(X) + CC 


LNiieGe. MUETIPLE.OF PI 


ARGP = AFG/PI 
Pins = TetereaRGP) -— 1 
IT = PABStITHTA) 


ae OR A NEGATIVE 


10 = MOD(IT.2) 

Wmwl@eF&O.0) 12 = -1 

THTA = PI¥DFLIAT(ITHTA) 

fpmNeseGi.2) THA = "TATA 90.5 Do2P 1 

HE VALUE OF X CORRESPONDING TO THIS INITIAL VALUE OF 
CALL XZERO(THT AsX1A+BsCC) 

NOG) eax 


mz 
Oo 
a 


A Lae PREG Ss Ohh > Kero UNITE SO8 ZEROES HAVE 


DO 2 J=2,500 
LHIA = Than + -P] 
XS hh SP eB) 
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CARE XZEROCTH Aa Xs Ayo CC) 

XOCJ) = X 

Deetove Gia ALh) GO TO 3 
2 CONTINUE 

J = 500 
3 CONTINUE 

RETURN 

END 


SUecoUl INE X22 RO (THIAy x AybeCC) 
Pees oOLVES THe EQUATION THETA = A®X + B*LN(X) + C FOR X BY 
AN ITERATION PROCESS 
Rie el Rees Cat Ag 7) 
eee eee CC St Oe) = 7 AA +B) 


Poet ree AR APPROMIA ION RES@LIS IN A NEGATIVE ESTIMATE 
ore Cte te tnOOmParts SOSREPLACE THe eslTIMATE Ur _x%. BY ONc— 
Dome © THE IMWPIITAL VALUE OF X AND REPEAT THE APPROXIMATION 
meee > A PeMe EST inate. REPEAT THIS PROCESS UNTIL THE 
SeresAl= CECOMES POSTTIVE 

Wet Wee 020002 XxX) = 0. 1s x 


Re eno oman ewe tee tO -—te) GO 10 2 
X = Xl 
CODToO™ 1 
2X = Xl 
RETURN 
END 
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SUBROUTINE CSINTCICSsMs,I1Q,X0eF oP) 


Mol? GAxRIES OUT TdHe Se ene PRP eGRAT LeN OF -7 He. BOUND- 
meee PRYTEGRALS FOR Z 1 


teem cl tT REALS S CA-h,O0-Z) 
EOrornal FPONPTPCN Sa STi ocos 
COMMON/MONE/C,0,FFA 
COMMON/ INT /XSM,FFs BUsXLM 
COMMON/GEZRO/5A,B, AMP, ICC, IQ] 
COMMON/RHO/N 
DIMENSION STS(4)-TI03),AC(4)+,AD04) ,AES(4) 5 XOCM) 
yec = ICs 
US ="ORCD0 
SINT = 0.0D0 
AES = Ona) 
DO 4 LL=1,3 

+ WRES(L E+) = ASS{LLY & P 
FN = DFLOAT(N) 


Teun GT sb Ss = PPelTEN = 21D) -UrN=+0.00)) 

AL = 0.25074FN + BO 

Priel eEs. 1) GO 10 7 

DO 5 J=1,4h 

ITFCXO(J).GE.XS) GO TO 6 
5 CONTINUE 
6edie=eles— 1 

LlsiteetsOveI) = 1 

ie Sine ER TheerU eT ION IS SWARTING A POSITIVE 
VE SWING ON @fHE NEMT HALF-CYCLE 
ITQI = MOD(JJ,2) 
DE Olen) 101 
IFC IT Ol2&O.0) IQ! 
Gb 10°36 


TQ 
=L¢ 


Isis 





100 


10 
5 


C AND 


ur 
es 


2 


Elon 


JJ = 1 

IQ] = 10 

D = XO(JJ) 
itis sO. 1 ) 
DO 9 LL=1,4 
Prt0.60.0.000) GO 10 100 
Aeree) = DEXxXP(—-AES (LE) =D) 7D 
GO TQ 9 
AD(LL) = 
CONT INUE 
DO 10 J=1.M 

oo) =f) Lg 

Pe UXON 1 SeT. AL) 
CONTINUE 

jal a ee er St mere ee ee ait 
DOM s 51 


oO ies 


Je0D9 


GO°TC ii 


Dey REPRESENT SJCEESSIVE ZeERn@eS OF THE TRUE 
= D 
P= -XOCK 41) 


INTEGRAND 


Tie Peeve eer Oo tie omer UNCTION WITH ZEROES AT C AND 


FA = 
FFA = 
Bio = 


Ome ai 7203 oo et o301/ (D> = 
FA 
DFLCOAT(1Q1) 


G) 


ane = Onleiee (0.100 —piIg 
IF(N-EQ.-1) GU TO 14 

DO 12 LL=1,4 

AC CUR) =) HOLL) 

Aoi) = DEXPL—AES (LLY *0I7D 


He iN Wee. Al JF THE APPROX BMATE FUNCTION IS EVALUATEC 
DO 13 LL=1+4 
Be > SRS!) 
REG = ACWW) 


SZ 





AED = AD(LL) 
CAEL CSSUR(D.+C,FA; TERM, AE,AEC, AED) 
Wer STS(LL) = TERM 


io =. See DOs we oS tl) — 0.301 (SiSt2) 
Ls MS Sg (Eo mes ae ee, 
GGy Te “16 


14 UL = Q.10D1 


ee APPROXIMATE FUNCTION MUST ALSO BE EVALUATED 


— 1 
OI 
> 


DO 15 LL=1.3 
BL = UL 
Cie = io bets P 
CALL DOQG32(BL+UL+FQN,TT) 
15 TI{LL) = TT 


is = ss = DIQ*AMP*XFA*(TI(1) - O-2D1*TI(2) + TI(3)) 
meee) «6©EVALOUATE THE INTEGRAL OF THE DIFFERENCE BETWEEN THE 
TRUE INTEGRAND AND THE APPROXIMATE ONE NUMERICALLY USING A 
32-POINT GAUSS-LEGEYDRE QUADRATURE 
ie OD = DFLGAT({] + (1 - IQOI1)/2) 
CC = peer o. 1) 1 
Pees thee AASeeGe [ter eeoOx IRATE FUNCTION, SO CHOSEN TO 
MAKE THE FUNCTION! BEHAVe PROPERLY 
B = FA=(0D"D - CCC) 
CAL B@GS2(C.0. FCN, FG) 
SPAT = SINT + FG 
17 IQI = -IQI 
Peel L Y THE REGION =-ROM J 10 XS AND FROM XL TO 100 IS INT~- 
FGRATED WITHOUT USING THE APPROXIMATE FUNCTION, UTILIZING A 
512-PCINT GAUSS-LEGENDRE QUADRATURE 
testes .6%.2) GO 10. 18 
oO = Cie) 


ii OsEomogomoy. GO 1024 
CALL 00512(0.0D0:D;DSIN3ZI) 
Corre 25 

24 ZI = 0.0D0 

25°C = Keer 1) 
i (Cer somes GG ib 0 
CALL DO512(C ,0.6D2DSIN, UI) 


13 





Ci = "0 602 


GUpepels 2 1 
ZomCLo= C 
UI = 0.ODO 


Peres OE «Oslo CO Tazo 

21 CALL DQ512(CL,0.1D3,;O0SIN;UJ) 
GO TO 19 

18 D = XO(JJ) 
IF(D~EQ.0.0D0) GO TO 28 
CALL B@512(0.)000,0,D0C05, 21) 

Comte =c 

28 Z1 = J.9DO 

27 C = XO) ) 
PE(C.GEsOmep2) GHeTO 22 
CALL DQ512 (Gse.GD2500C0S,U1) 


wr = 0.602 

@O 10 23 
22 CL =C 

UI = 0.ODO 


PriGeGe-.Os1 02) GO TO 26 
22> MALL DOQS120CL,0.103,0C005,UJ) 
eo TO i 
26 UJ = 0.0D0 
ges INGEGRAL ES THE SUM OF ALL THe INDIVIDUAL INTEG- 
RATIONS 
Lo FOS Svat cl eee) ee 
RETURN 
END 
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7 


C 
B 


S 
E 


5 
5 


UM FINDS THE INTEGRAL OF 


nec 


lJ 


SUBROUTINE CSS UM(D,CsA,TERM,AE,AEC,AED) 


(ah) > 


ERA oo CCE SoOrVe cero o” Or tie tN 
PAPC CT Sees CAS eZ) 
CQHMHMON/RHO/N 


LIM 
DEN 


S = 
Fo= 
DI 
DR 
PHP 
FN 
AC 
AD 


oo 
— 


Cd 


Nest 1 
~~“ DSO rime A eae 7 
DATAN2(A,4AE) 


Ie DDI 


0.101 


O.1D1 

OWL DE/DE 
PH=DFLCA T(N+2) 
Ose Dle/ BEWOAT dale Wh ) 
AEG 

AED 


— 


— 
-_ 


DOF UGA TC) 

PN F I 
DI=DEN 
be eDEM 


EXE] 


PHP — P4 
AC*C 


(E 
TE 


IF( (AC. EQ.0.909) -AND.(1.EQ.1)) 


TERM 


S 


AD*D 


Ale 
EGR 


Le 


(-ALPHA=X))*SIN(AX+B) 
RAWND 


= 0.101 


(ee Se CRT laD Rat AD + AC)/F 


= Files /.D1 


RETURN 


END 
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SUBROUTINE COGSZ(XLsxU,FCT.Y)} 


DOG32 IS A 32-PDINT GAUSS-LEGENDRE QUADRATURE WHICH INTEG- 
Pat ES “\)  liiiaaaias Ceara tee Ge.  SEIWeEEN THE LIMITS 
XL AND X 
DESCRIPTION OF PARAMETERS 
XL: LOWER BOUND OF THE INTERVAL 
ROU PSR we OUD: CF JHE INTERVAL 
Pel. aie OfeetHe we xXPeERtALLY DEFINED FUNCTION 
TI>BEs (INTEGRATED 
Y: RESULTING INTEGRAL VALUE 
DOWBLE PRECISION Xb sXe A: By GarGGl 
Der TNTMG A AND B ENABLES THE VARIABLE TO BE CHANGED TJ CON- 
son Wate let be tea Oe ool Eee eur THE GAUSS—L EGENDRE 
U 4 ‘es 


A=. 500% ( XU+XL) 

B=XU-XL 

C=.498631930924740/8D0*B 

Voeo nee e041 Soe ese e Cla te CT CA=C) ) 
C=ee@Sa@e0575577 2654 1.70 0-6 

Yoweri si lotsa coon ee (FCT IAPC} +FCTCA-C)) 
C=s4Oesoli ec I ot ooe2eO 6 

Y=Yt. 1 26903205463 103500-1(FCTCASCI+EFCT(A-C) ) 
C=.46/45303/795886984D0B 

Your. 1/1367 oReeeo10(17/0-1* (EC are aeECitA-C)) 
C=.4481635/7883)26 J6D0=B 

Yo vite d tne eee 1) 3400-1 Clee tea CAC) 
C=.42468380 686 62 8499) 0*B 

YHYt. 2547332505) oem oU-) (EF Cimbate AEC? (A=C) ) 
C= 6.391 24189 (96355 /12600=8 

V1 eco oc O tee 1 (PCY CAG) +eCr (A—C) ) 
C= .3650910593/014484D0*B 

Yaar eto COL BOSS (ECT tUATC)+eEC1{ A-C)) 
C=. 3315221334651 976900 *B 

YH¥+. 3617289 1054424253D-1*(FCT(A+C)+FCT(A-C)) 
C=. 295 51618020501 16098 

NS re 0541592550 lS BUH 1S (Pt (AC) + ECT (A=C ) ) 
C= .25344995446611470D0*B 

Y=¥t.4165559621] 13473378D-14 (FCI (A+C)+FCTI (A-C) ) 
C=.2106756280653176700*B 





PON EV AlWATES THE AP? PROXIMATE FUNCTION AT GIVEN 


Yoy+ 3 826 046 Dee 2m SIGD-1¥(FCT(AtC)+FCT(A-C)) 
C=.16593430114106382D0*B 

ett s 4 oS Cos 39547 eer 2420-1 * (PC T( ATC) +FCT(A-C)) 
C=.119643681125)6854D0*B 

¥ = ¥ eos 22 199540402 2830-1*(FCT(AtC)+FCT(A-C)} 
C=e(eeeDpI 80791390 825D-1+5 

oy eto los oO) S763 1 4o0D—- i= tr Cli Ate y+ECiCA-C)) 
C=.241 53 Sav 843 869158D—178 

VaR YP. 460217394425 7363900D—1 mr FC! (AFC )+ECTCA-C))) 
RETURN 

END 


DOUBLE -PRECTSTUN Pagel LON=roNtrZ) 


We liciyYy ReAc* 8 (A-H,0-Z) 

COMMON MCNE/C, 0, FA 

EChe=-tDEwet-o=2) + DSKP(—0e72))7(2-2 + FATFA) 
RETURN 

END 


DGUSLE “PRECIO Forme Ten aC NieZ) 
“GIVEN 

We ec eee 8 6 (A-f, 8-2) 

CUMIN? CONS'7 PHl,.PH2,F RV PFRZ,ASCL,P1,CA,CI1,AA 
COMM INYGEZRO/-4,86;,F,1CS,.1 0] 

COMMON/ RHO/N 

COP OMe eee 2), ICS 


FR = FR1 
PH = PH1 
Peck 
GO TO 3 


Cee =) FiZ 


157, 


(digg) 101) eG | Orca k 


NCE eso THE TRUE INTEGRAND AND 





Fg O's 2D. 

ARGF = FR*¥Z + AA*DLOG(Z) + PH 
PeaGH = FPA=Z +98 

See OU. 1B CSP - he) 

a= Oe oly (22) 


NN = N + 2 
DO 6 I1=1,NN 
ZN = ZN*Z 


PNV = S*S*S=DEXP(—Z)*ZN 

PF (1CS.GTw2) GO TO 4 

FON = ENV*(DSIN(CARGF) ~— F*DSINCARGH) ) 
oo TOS 

Pen = ENVE(DCOSCARGF) -—- F=OSINCARCH) ) 
RETURN 

END 
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SUBROUTINE DQ3512(C,D,FCT,Y) 


DP ieeClT REAL*8 (A-H,O-2Z)} 
COARON/L51 2/4256) 741256} 
COMMIBNY7 GON ST/PH1+PHZsFRI1sFRZ2,ASCL,PI ,CA,CT;AA 
COMMON/GEZRO/FA,B,AMP, ICS; IQ] 
COMMON /RHO/N 
GO. TOS ls24 ic ) yl Os 

nets rx FR1L 
Bro = Pid 
P= SASL 
GOsTO 3 

Zee R. =ark2 
Pr = sPad2 
P = 0.1D1 

3 Y = J.9D9 


i THe VARI ASEE SO THAT THE INTEGRATIGN LCIRITS BECOME 


DPC = 0.50D0"(D + C) 
pee = 0.50080D - C) 
DO 4 1=1.256 
XDC = Xt 12 ore 
ZP DEC. taexDe 
Life = OPC = KDC 
ZMN = O/1D1/(ZM*2Z) 
ZPN = 0.1D1/(ZP*7ZP) 
NN = N + 2 
DC 5 J=l,NN 
ZMN = ZMN*ZM 
5 ZPN = ZPN*XZP 
SP = 0.1D1 - DEXP(-P*ZP) 
SM >= 0. 1D)L.— MEXP(—Pmz Mo) 


Eye 





SMe = Sr sheer DEXR{=7P 4 2P ly 

ENV = SM*¥SM*SM*=DE XP (-ZM)*ZMN 

BEGP = FREZP=+ FR=OLCG UZP lee PH 

ARGM@ = FR¥Zhi + AA=DLOG(ZM) + PH 

¥ Y Sean Te -CENVP FC LUARGe) 2 ENVN=FCT CARGM) ) 
y Yee 

RETURN 

END 
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SUBROUTINE Hi Vrs! ;C2yo55CC , I MAX MAX } 


MeeisP EVALUATES THE BOUND-FREE INTEGRALS FOR THE CASc Z=1 


1 


PAPCI CIT REA eA Hse: 
EXTERNAL FACT 
COMMON/CONST/ PH1,PH2,FR1,FR2,ASCLyPI,CA,CIsAAA 


ase RHO1(4),RHO2(4),THTAS1(4),THTAC1 (4); 


THT4S2 (4),THTSC2 (4) 
DIMENSION S1(IMAX),C2¢( IMAX),SS(MAXI),CCCMAXTI) 
Fl = FRl**2 
EZ =, Phe ae 
DO 1 K=1,4 
A = OF LOBT(K - 1) 
RHO1(K) DSGr Tl OO. wbie +A SC 2 et EF 1) 
RHO2(K) Deen Titel] + Alt=2° 4 FZ) 
THTAS1 (K) Gers IN CPR REO TK) 
THTACI1(K) Oar OSI (0, ID] + A=ASCL) /RAG1L CK) } 


I 


ThtAS2t0 KR) = OARS INGCERZ/RAOZ CK) ) 

THEACZ(K) = CART OStCO.IDI + AY 7RHEZUK) ) 

Sl(1l) = Cee ee eo wy (0..101 + ASCL)) 
1 — DAT Mien 1/7 (OsTO] + Os20i7-ASCI))) 

2 = DATAMSRI7 tO. IDL + Jo 0)AS5Cu}) 

C201) = QO. opie UOC Wet Oo 1D et Ase) 2er.F1)/( (0.101 
1 + On 2 OS Clyrrne ot EI) )) tee oD DLOG((( 9.101 
2 + (OsSDieea sc) ae + Fly 7 COS + 1) ) 

SS) = DATA R2) 0.3 Dias AT Te ERZ/9~.201) 

Js DATAN(FR2/0.301)) — DATAN(FR2/0.4D1)} 

CC(l) = Unsere LUGO. 40) + F2)7 Vaso e F2)) + 
i Ore PomaLOGt(O.16GD2 + F2Z)/7/1(0.1D1 + F2)} 

DO 2° 3225 

= J soa 

Fl o> VELGAwel) 

ECT = Fae. 2) 

S143): =e OSI NCE I= ties) AeA) =I] 

i = 07 SDT UD SN it Ser yee eta jens | 
2 = Use bis) ARO bt Ss) =] } 
S SO Sere rome) Rt) Olean 1 ) 

C23) =SFGi= (UDECOs nites e (hy RHO Cl) =] 

l =) Oe 3D T=40COSt( EI =IBIAC II 2))/ RHO] ( 2) 1 
a = TOC SE tenes) 17 RAO 3 eT.) 

3 = WCUSURIIRIACI (4) )/ RHO) (4 ) =e! ) 


Oe ero a) We) a Ce 
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2 CONTIG E 


RETURN 
END 
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SUBROUTINE A1L2SUM(NUsLMDA+AAVAL) 
EVALUATES AL(NUsLAMDA+2) OR A2(NU,;LAMDA?2) USING I] 
Reoreepivecy 

IMPLICIT REAL®8 (A-H+O-2) 

DIMENSION VAL(1) 

COMMON/TRM/ALFA?ETA 

COMMON/ATNO/Z 

MLFA2 = Z + ALFA 

FCTR = 0-.2D1L*ALFA2/(0.2D1*Z + ALFA) 

Ww) = NU + 1 

he = NU + A 

ee = COMPO ORROAT (NU 

ECt = 2O0esbiy ECT 

TVAL = O0.ODO 

DO 10 JT=1.NU 

Fo = PRDEUBMAT(IVNUL - T) 

POT = PO rere 

TVAL = TVAL + FCT=VALONL-I)*F 

AA = TVAL/ALF42¥%*(NL-2) 

RETURN 

ERD 


sie che Be oe se ote oie che ois ok ok ok fe ote oe ate oe oe ee oe oe ok eo OE MC SE OE Ss Ae EE OK sie deck se ke ds ake fc ok Fe OK OK At 
se ese oe ak ok ak oe Sk a ate eC RCC SK a BC ia EB Bea a a CH AR ae a ie ee a ee IE BK EZ «SK US OK IG 
we dec we ae hg see se eS oe Se STK Be Ha ak ek ca ok aK ae ake ea yak a ok ak ak ak aaa ak sa a 
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SUS @UTINE PHASE(AS;AC,8S,BC,HAXI, HS*HCyE) 
pemeotri E@Eite fom Sr E)VCHIt)y> AND 


Deel Reabees (A—-H,y0=Z) 
COMMON/TRM/ALFA,ETA 

COMMON/ATNO/ Z 

COMMON INDEX/NIYNJ,LIsLJeMIsMJ,IS 
DIMENSTON AS(1),AC(1),BS(1),BC(1) 

DIMENS TON NU(26),LMDA( 26) ,MU( 26) ,COEF( 26) 


FN = DFLOAT(NS) 

FL = DFLOAMEL ) 

Eie= DFLCA Limi) 

peeveUe sO anes, LANUAg c UR AND Tie SCOEFFICIENTS FOR 
Pe 26 Siew OT Sie era ih in EWEMENT 

ee) = beet | 


Nut 3) = NUCL) 
NU( 4) = NUCL) 
NUC 6) = NUCL) 
NUt &) = NUCL) 
NUC 9) = NUCL) 
NU(19) = NUCL) 
LMDA(12) = NUCL) 
LMDA(26) = NUCL) 
NU{ 2) = LJ 
LMOAC17) = NUC2) 
NU(21) = NU(2) 
Mi) = ot 2 
NU(20) = NU(5) 


NU(22) = NUC5) 

LMDA(10) = NU(5) 
LMDA(11) = NU(5) 
LMDA(14) = NUC5) 
CMDAtLS) =e 5} 
LMDA(16) = NU( 5) 
LMDA(18) = NU( 5) 
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rma (24 ) 
NUC 7) = 
NU{10) = 
NU(1l2) = 
NU(13) = 
NU(15) = 
NU( 17) = 
NU(18) = 
we( 23) = 
WMDA{ 3) 
Cea 2 2 ) 
NUC11) = 
WW 25) = 
LMDA( 8) 
NU(14) = 
NU(2Z4) = 
NU(26) 
LHDAC 1) 
LMDA{ 2) 
EMeA{ 5) 
LMDA( 6) 
LMDA( 7) 
LMDA( 9) 
LMDA( 20) 
NUC16) = 
LMOAC 4) 
Liter 19 ) 
LMDA(21) 
EMmDA(1 3 ) 
LDR 23) 
LIMDA(25) 
MUC 1) = 
‘Ut 2) = 
MU({ 3) = 
MUC 7) = 
MU( 8) = 


= NUC5) 
ee: 1. 
Nei «] 
NU( 10 ) 
NU( 10) 
NU( 10) 
NU( 10) 
N12) 
NUC 10) 

= NU( 10) 
= NUC10) 
NJ 
NU(11) 

= NUC11) 
Nek. 2 
NU(14) 
NU 14) 

= NUC14 ) 
= NUC 14) 
= NUC14) 
= NU(14) 
= NU( 14) 
= WUC 14) 
= fU(14) 
Nom 1 

= NJ + 3 


= LMDA(4) 
= LMIA(4) 


= bot 3 


= VMIA( 13) 
= LMBA13 ) 


ORE Ee 
MU(1) 
MUC1) 
MU( 1) 
MU (1) 
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MU( 10) 
MUC 11) 
MOC 12) 
MU(16) 
MUC17) 
MU( 4} 
rel eee) 
MU( 6) 
MU( 13) 
MU(14) 
M15) 
MU(19) 
Mu t2Z0) 
mu C21) 
MU (22) 
MU (23) 
MU024) 
mocz 5 } 
MU(26) 
MU( 9) 
MU(18) 
COBF( 


COsrt 10.) 


Comer 


CULM REZ) 


CGE 


COEPOI)) 
CGEE lie.) 


COEF ( 


COEF (13) 
COER CES) 
COEF(19) 
COEEtZ 0) 
CGiErIT 23) 
COER (24) 


COEF ( 


MU(1)° 
tO 
MU(1) 
MUC1) 
MU (1) 
J 
MU(4) 
MU(4) 
MU(4) 
MU (4) 
MU(4) 
MU(4) 
MU(4) 
MU{ 4) 
MU (4) 
MU(4) 
MUC4) 
MUC4) 
HUC4) 
MJ + 1 
MU(9) 


1 


i 


i 


=9. 25D0@WLF Aawee — E 
COEF ( 
O.5DO¥ALFA¥(FL + 0.101) - 
COB 
O.5D0*ALFAX(FN 4 
Geet 3) 

0.5 DOXALFEA“FM 
COEF (4) 

COert 44 

COEF(4) 
Se clee tt 
Sees ten 
Sweet 
Omen) 
—FM% (EM 


1 


ag) 
Cae) = 


4) 


a eh te alee ge ne Qe OT 


166 


Z 








COEP. 15) =-COEF( 6) 


CHEE. 7 )a=—-0s5002FE (EL On1U 1} 
COEEtl 7) = COSE(= 4 
COEr( 8) .= -Os 5002 Fir 2 Dy) 


COE GIG) = COEPC Re) 
COEFY 9) = 0. lybd 
COSmyrc) = CORT >) 
COERpO2Z1) = feet 
COERCZO)@=6COERIZ1) 
COERG22). Sebi 
Cire 25)7 =. COEPtCZZ) 
We( 1S .—EO.0) 63.710 16 


maeGe THE STGNS OF TRE UASP is IERWieerGR THE TRIPLET CASE 
DO 1 I=10;18 


L“COEPUL) = =COErG 
GO) 2al= 2 326 
. COEPII) = -CUOER a 


LO wEaaAK = MAN! 
NNMAX = MAXI + J 
HS = 0.000 
HC = 0.0D0 
DO 20 K=1,18 
IN NUK) + 320+ ais Xe eA) + GM AX=MUCK) - 1) 
HS = BS 4 See) COEr CK) 
eC = oie + GC Ce) EO BE) 
20 CONTINUE 
SO SUM 9 426 
PNP =NUTK) Fo 2 + NPA X= COMOAVK) + LLMAX*MUCK) - 1) 
pS =H So+ BSCTN) COEF (K) 
HC = Gere eC (ON iC OE EAGK 
30 CONTINUE 
RETURN 
END 
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13. 
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>. 
EG. 


ves 
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